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SUMMARY 


Approximate nonlinear inviscid theoretical techniques for predicting 
aerodynamic characteristics and surface pressures for relatively slender 
vehicles at supersonic and moderate hypersonic speeds were developed. 

Emphasis was placed on approaches that would be responsive to conceptual 
design level of effort. Second order small disturbance and full potential 
theory was utilized to meet this objective. 

Numerical codes were developed for relatively general three-dimensional 
geometries to evaluate the capability of the approximate equations of motion 
considered. This report represents a user's manual for the second order 
analysis and optimization codes. Contained herein for each program is a brief 
description, its input data, variables, and subprograms, a flow diagram, and a 
sample test case. Results from the computations indicate good agreement with 
experimental results for a variety of wing, body, and wing-body shapes. Case 
computation times of a minute on a CDC 176 were achieved for practical 
aircraft arrangements. 
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INTRODUCTION 


An examination of the literature for supersonic/hypersonic aircraft 
provides an indication of the flexibility and generality required for an 
analysis technique. Typical configuration development variables include wing 
section, incidence, height, dihedral, planform, effectiveness of longitudinal 
control surfaces for trim, effectiveness of empennage for directional 
stability, and propulsion system-airframe interactions. 

Response to these requirements at the conceptual design level have been 
by the hypersonic impact methods and the linearized analysis and design 
algorithms. These approaches can treat complex geometries efficiently with 
minimum response time and cost. Shortcomings exist with both methods. The 
impact methods ignore component interference effects and crossflow 
interactions. The linearized approaches exclude shocks, vorticity and entropy 
wakes and layers. 

' A need exists for an improved analysis technique to optimize vehicles 
designed to travel at supersonic/hypersonic speeds. The technique should be 
more accurate than simple noninterfering panel methods and more 
computationally efficient than the current explicit finite-difference methods. 
Enough of the physics of the flow should be included to allow realistic 
optimization and permit considerations of component interaction. Nonlinear 
potential theoretical formulations hold the promise of meeting these 
objectives for preliminary vehicle definition efforts. 

To satisfy the analysis needs, a program was undertaken to advance the 
aerodynamic prediction capabilities at the conceptual design level. Numerical 
second-order potential small disturbance analysis was developed as a first 
step up from the widely used linear theory. Such formulation incorporates 
nonlinear behavior by iteration of the Prandtl-Glauert approximation. This 
approach is known to extend the prediction success for airfoil and cone 
surface pressures to substantially higher values of the hypersonic similarity 
parameter than the first-order theory. References 3 and 4 contain the details 
of the theoretical development of the second order analysis and optimization 
method. 

This report provides a user’s manual for the second order analysis and 
optimization codes. Contained herein for each program is a brief description, 
its input data format, a complete list of variables and subprograms, flow 
diagrams, and sample test cases. 
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The computer program entitled "SOPA-Second Order Potential Analysis and 
Optimization Program" can be obtained for a fee from: 

Computer Software Management and 
Information Center (COSMIC) 

112 Barrow Hall 
University of Georgia 
Athens, GA 30602 
(404) 542-3265 

Request the program packgage by the designation LAR 13314. The programs are 
written in FORTRAN V for use on the Control Data 6600 and the Cyber series of 
computers. 
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ANALYSIS PROGRAM WBODY 


DISCUSSION 

Computer program WBODY performs a linearized analysis of wing body 
combinations using axis and surface singularities. The solution satisfies the 
Prandtl-Glauert equation with boundary conditions as prescribed by the 
assumptions of thin wing theory. A second order solution may then be 
performed using the results of the first order (linearized) solution. The 
second order solution on the lifting surfaces is performed by using a 
modification of the exact result available for flow in two dimensions. On the 
body the second order solution combines an exact axisymmetric axial solution 
with the first order cross flow solution. The second order solution is not 
valid for transonic Mach numbers and should not be used for Mach numbers 
between 0.70 and 1.60. All output is placed in a dataset which may be used 
for computer graphics. 

If a body is present and the use of axis singularities is indicated, 
isolated body axisymmetric axial and cross flow solutions will be performed 
first, using only the axis singularities. The axis singularities consist of 
linearly or quadratically varying line sources and quadratically or cubicly 
varying line doublets. The singularity and control point spacing are input by 
the program user. Either mass flux or velocity boundary conditions may be 
specified. A second order axial line source solution will be performed if 
indicated by the input. 

The body panels are quadrilaterals of arbitrary shape having a constant 
source distribution. The lifting surface panels are quadrilaterals having two 
streamwise edges with a constant vortex distribution to simulate lift, and a 
constant or chordwise and spanwise linearly varying source distribution to 
simulate thickness. 

The first order solution is performed by satisfying the condition of zero 
normal velocity, or zero mass flux if desired, at the panel control points of 
the body and lifting surfaces. The lifting surface boundary conditions are 
linearized and are satisfied on the mean camber line. The normal velocity or 
normal mass flux is composed of the sum of contributions from the free stream, 
angle of attack, and perturbations from the axis singularities, body panels, 
and lifting surface panels. 
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The complete first order solution is composed of a linear combination of 
four basic solutions: 

1. The cross flow or add load solution. This solution contains 

no lifting surface thickness or camber. 

2. The camber solution. This solution satisfies the boundary 

conditions for the lifting surface camber with zero free 
stream velocity and zero normal velocity (or mass flux) 
on the body panels. 

3. The lift due to thickness solution. This solution represents 

the additional load due to the normal velocities induced 
by the thickness of the lifting surfaces. 

4. The lift due to the body solution. This solution represents 

the additional load due to the normal velocities induced 
by the body in axial flow with no camber. 


The second order solution, if desired, is performed using the results of 
the first order (linear) solution. 


INPOT DATA 

Data is input using subroutine DECRD1 , described on page 19, and is 
stored in the array called "DATA". All locations are initially set 
equal to 0. Each case contains: 


Miscellaneous data locations 
Body geometry data locations 
Lifting surface geometry locations 


6-88 
701 - 1900 

1-5, 89 - 700 

and, 1201 - 1900 


All of the data except the lifting surface geometry data must be input first, 
with the last card having a negative location number. This data must include 
any input for body panels and axial singularities. 

The lifting surface geometry data follows with the last card for each lifting 
surface having a negative location number. 

The last surface or last card for a given case is indicated by a positive 
value in location 1. 
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The body or surface geometry may be read from an APAS (Aerodynamic Preliminary 
Analysis System, references 5 and 6) file using data location 14. 

The first card in each case may be a title card containing up to 72 
characters. If the first card contains a blank or a minus sign in the first 
column and blanks in columns 2 thru 4, it is assumed the card contains 
numerical data, and is not a title card. 

In addition the first card for each lifting surface may be a title card 
containing up to 16 characters, i.e. cols 1 - 16. This title will be used 
for the name of the surface. 


The following is a description of the possible input for each case. 

LOC VARIABLE DESCRIPTION 

1 <0. The run will terminate. 

> 0. Indicates the last surface has been read. 

> 1.1 The run will terminate after the slender body 

geometry (if any) has been printed. 

2 N The number of panels in the chordwise direction for 

this surface. 

3 M The number of panels in the spanwise direction for 

this surface. 

7 ACALC 2. Compute the aerodynamic influence matricies and 

the quasi-inverse matrix and store on unit 11. 

1. Compute source influence matrix and store on unit 
11. The vortex and quasi-inverse matricies are 
assumed to already exist. 

0. Compute new aerodynamic influence matrix. 

-1. Use aero influence matrix stored on unit 11. 

-2. Read aero and quasi-inverse matricies from unit 11 

(A second order solution always requires a 
quasi-inverse matrix). 

8 x.y If x.y < 0. a 2nd order solution is performed. 

If y > 6 d/dy of v terms are included. 
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ITYPE the type of source panels used. 

If =-2. linearly varying source panels are used 
(spanwise and chordwise without leading 
or trailing edge panels). 

If =-l. linearly varying source panels are used 
(chordwise varying only without leading 
or trailing edge panels). 

If - 0. Constant source panels are used 

If = 1. linearly varying source panels are used 

(chordwise varying only with leading 
and trailing edge panels). 

If = 2. linearly varying source panels are used 
(spanwise and chordwise with leading 
and trailing edge panels). 

Only the geometry and the axial 

singularity solution will be calculated. 
The program will stop after all of the 
geometry (including axial singularity 
geometry) is calculated. 

13 CENT 1. Control point at panel centroid. 

Control point at panel center otherwise. 

14 <1.0 Obtain body and aero surface geometry from input 

1.0 Obtain body geometry from APAS (panels) 

1.1 Obtain slender body geometry from APAS, 

create panels, and save after aero geometry. 

1.2 Obtain slender body geometry from APAS, modify, 

create panels, and save after aero geometry. 

2.0 Same as 1.0, 1.1, and, 1.2, but includes 

2.1 aero surface panels and geometry. 

2.2 

2.0 Use twist and cambers from APAS and input data. 

2.1 Use twist and cambers from input data only. 

2.2 Use twist and cambers from APAS dataset only. 


^ 3.0 Obtain only aero surface panels and geometry from 
APAS (no body). 

3.0 Use twist and cambers from APAS and input data. 

3.1 Use twist and cambers from input data only. 

3.2 Use twist and cambers from APAS dataset only. 


12 If < 0. 

If < - 2. 
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IPRNT 


O.K 


INTMED 


DPLOT 


XIJ.XKL 


XIJ.XKL 


1. Print source dz/dx and z/c matricies 

2. Print both source and camber arrays 

3. Print camber slope matrix 

0. No printing of panel geometry. 

1. Print body panel geometry (source panels). 

2. Print vortex panel geometry 

3. Print body panel and vortex panel geometry 

-1. Print body panel geometry from APAS dataset. 

-2. Print vortex panel geometry from APAS dataset. 

-3. Print body panel and vortex panel geometry (APAS) 

K > 0 Print vn due to thickness. 

K > 1 Print u due to thickness. 

K > 2 Print v due to thickness. 

K > 3 Print w due to thickness. 

K > 4 Print phi due to thickness. 

Various orders of intermediate printout (-1. to 4.) 

-1. Least printout (no upper and lower surface Op's) 

2. More Printout 

Prints add load solution and v-normal 
(v-w), (w-b), (b-w), (b-b) 

3 . Above + 

4. Above + 2nd order boundary condition solutions 

5. Above + odd and even 2nd order velocities and Cp' 

0. No data is placed in a plot dataset. 

1. Data is placed in an APAS dataset (geometry) 

2. Data is placed in a DDP dataset (for plotting). 

' 2. Data is placed in a UDP dataset (extended) 

DPLOT > 2. is required if the dataset is to be used for 
a 2nd order optimization calculation. 

Detailed influence matrix calculation printout for 
vortex panels. From subroutine VORTEX. 

XI J * 3 digit number for control point 

XKL = 3 digit number for panel number 

Detailed influence matrix calculation printout for 
source (thickness) panels. From subroutine PHIS11. 

XI J = 3 digit number for control point 

XKL = 3 digit number for panel number 
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APRNT 


IJ.RKR 


.ne. 0. Print the aero influence matrix (vortex and body) 
.gt. 0. Only normal velocities are printed. 

.It. 0. All velocities and phi are printed. 


I 

or J = 0 

nothing 


I 

or J = 1 

vortex 

I - influenced 

I 

or J = 2 

body 

J - influencing 

I 

or J = 3 

both 



11. RRR 

vortex on 

vortex only 


12. RRR 

body on 

vortex only 


21. RRR 

vortex on 

body only 

ABS( APRNT) = 

22. RRR 

body on 

body only 


13. RRR 

all on 

vortex 


31. RRR 

vortex on 

all 


23. RRR 

all on 

body 


32. RRR 

body on 

all 


32. RRR 

all on 

all 


The influence matrix printout will stop after RRR rows are 
Printed. If RRR = 000 all of the rows are printed. 

24 SPRNT * IJ.RRR 


.ne. 0. Print the source (thickness) influence matrix 
> 0. Only normal velocities are printed. 

< 0. All velocities and phi are printed. 


ABS(APRNT) = 


1=0 nothing 

1=1 vortex I - influenced panel 

I = 2 body 

I =3 both 

IJ.RRR source on vortex only 

2J.RRR source on body only 

3J.RRR source on all 


The influence matrix printout will stop after RRR rows are 
printed. If RRR = 000 all of the rows are printed. 
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25 > -1. Print perturbation velocities from axial singularities 

=0. Print velocities from sources and doublets 

= 1. Print velocities from sources 

= 3. Print velocities from sources (1st and 2nd) and doublets 

= 4. Print velocities from sources (1st and 2nd) 

> 4. Print velocities from sources (1st and 2nd) and doublets 

26-28 S.K Surface to be extended to centerline. 

If > 0 the surface will only be extended when the aero 
influence matrix is calculated. 

S * surface number (in ascending order). 

K * 0 Extend inboard panels with zero sweep. 

K * 1 Extend inboard panels with same sweep. 

K * 2 Extend inboard panels with negative sweep. 

i.e. if DATA(26) « 2.0 the calculation of the aero 
influence matrix for surface # 2 will be 
performed by extending the inboard panels to 
y ■ 0. with zero sweep. 

if < 0 the actual panel will be extended and locations 711-714 
should be used. 

29 < 0. beta * y / x is printed for each panel 

30 MN Mach number used for order 2 solution CR2 array. 

MN > 0. use MN as Mach number. 

MN = 0. use freestream Mach number (generally used). 

MN < 0. use normal Mach number. 

if < 0. the maximum value of CK2 = - MN 

31 X00 the x value of the pivot point for angle of attack. 

used for second order theory only. 

32 XCG For computing Cm, Cr, Cw, Cy 

33 YCG For non-symmetric rolling moment only 

34 ZCG For non-symmetric rolling moment only 
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35 MACE The Mach number 

36 RD The x/c fraction of each panel where the normal 

velocity is interpolated to calculate the zero 
suction drag. If 0. a value of 0.515 is used. 

37 RS The x/c fraction of each panel which is used to 

curvefit (Cp,x) in order to obtain an approximate 
value for the leading edge suction, (default = 0.250) 

Drag Polar (41 points calculated) 


38 

39 


Delta CL for drag polar. Default ■ .05 

Starting CL for drag polar. Default =0. 

40 

CBAR 

The reference chord length for computing Cm 
if 0. CBAR = CAVG is used. 

41 

CAVG 

The reference chord length. If 0. 
CAVG = SREF/SPAN 

42 

SREF 

The reference area, if 0. 
SREF = total area 

43 

SPAN 

Span, any consistent set of units may be used. This 
value is used for reference purposes only. If 0. 
SPAN = 2. * Y-max 

44 

RATIO 

The chordwise control point location, if 0. default is 
0.875 for Mach >1. 


0.875 for Mach < 1. 
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45 


Data locations 45 - 47 apply to 2nd order solutions only. 

IJK Printing of odd and even symmetry u,v,w velocities, 
if DATA (45) = LMN. I = L, J=M, K = N 


I J K 


* 0 

No 

velocities printed 

= 0 

none 

= 0 

none 

= 1 

u 

velocities printed 

= 1 

odd 

* 1 

add load 

= 2 

V 

velocities printed 

= 2 

even 

= 2 

twist and camber 

> 2 

All 

velocities printed 

> 2 

all 

= 3 
= 4 
> 4 

thickness 

axial 


e.g. To print odd and even symmetry thickness v velocities, 
use DATA (45) * 231. 1*2, J=3, K * 1 


46 IJK Printing of odd and even symmetry d/dz of u,v,w 

47 IJK Printing of odd and even symmetry d/dy of v 

I >0 is required 

48 J Printing of d/dx of camber and thickness normal 

velocities. 

J = 1 Camber only. 

J = 2 Thickness only. 

J > 2 Both 
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50 ALPHAO The angle of attack of the surfaces with respect to 

the x,y plane (degrees). Used for second order 
theory only. ALPHAD is made up of ALPHAO and 
the angle of attack of the freestream with respect 
to the x-axis, ALPHAI. 

i.e. ALPHAI = ALPHAD - ALPHAO 

For a first order solution only ALPHAD matters. 

51 ALPHAD(l) The angle of attack of the configuration with respect 

to the freestream (degrees). 

52 ALPHAD ( 2 ) 

58 ALPHAD(8) (maximum number of angles of attack = 8) 

If ALPHAD(K) > 90. an add load solution is perfomed. 

i.e. (ALPHAD = 1.0) - (ALPHAD = 0.0) 

61 CAMTHK(l) Input in form OB-OV.CAM-THR, and used with ALPHAD ( 1 ) 

OB * The order of the Cp on the body (one digit). 

OV * The order of the Cp on the lifting surfaces) 

CAM > 1 Camber included. 

THK > 1 Thickness included. 

CAM ■ 0 No Camber included. 

THK * 0 No thickness included. 

0B, OV, CAM, THK are each one digit. 

If ■ 0. OB - 4, is used for the body Cp, 

OV, is determined by DATA(8), and 

camber and thickness are included 

OV » 1, unless DATA(8) < 0. 

If OV > 2, OV = 2 will be used, and also the 

difference between the 1st and 2nd 
order delta-Cp's will be printed. 

If < 0. Same as above except the 2nd order axial 
solution (if performed) is used to 
calculate u,v,w. 
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The body pressure formula is 
determined by the value of OB. 

OB Cp 

1 Cp = - 2 * u 

2 Cp = - 2 * u - beta2 * u*u - v*v - w*w 

3 Cp = - 2 * u - u*u - v*v - w*w 

4 Cp ■ Isentropic pressure formula 

5 Cp = Isentropic for alpha * 0. + isentropic add load 

u,v,v use 1st order axial contributions if CAMTHK > 0 
u,v,w use 2nd order axial contributions if CAMTHK < 0 


e.g. CAMTHK = - 31.02 


indicates, 2nd order u,v,w from axial solution (if performed). 
3 indicates, pressure formula #3 on body. 

1 indicates, 1st order Cp on lifting surfaces. 

0 indicates, camber is not included. 

2 indicates, thickness is included. 


62 CAMTHK(2) used with ALPHAD(2) . 
68 CAMTHK(8) used with ALPHAD(8) 


The following data (89-700) are read for each lifting surface. 

89 TC The t/c due to thickness for this lifting surface 

If < 0. a thickness. influence matrix is calculated 
although TC = 0 is used. A second order solution 
will always calculate a thickness influence matrix. 

Locations 90-94 are the thickness coefficients 


90 

91 

92 

93 

94 


If 90-94 are all 0. a NACA 4 digit airfoil is used. 
If CO < 0. a biconvex airfoil is used. 


CO 

Cl 

C2 

C3 

C4 


The coefficient 
The coefficient 
The coefficient 
The coefficient 
The coefficient 


of the SQRT(x) term for 
of the x term for 
of the x**2 term for 
of the x**3 term for 
of the x**4 term for 


thickness 

thickness 

thickness 

thickness 

thickness 
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101 

SWEEPL 

102 

SWEEPT 

103 

ROOT 

104' 

ROOTXO 

105 

ROOTZO 

106 

SPAN0 

107 

XSP 

108 

FLAP 

109 

DIHDRL 

110 

RATIOY 

111 

SYM 

112 

ND 

113 

MD 


121 

Z(J) 

161 

Y(J) 

201 

YC( J) 

241 

XLE(J 


281 ZTE(J) 


The leading edge sweep in degrees. This value is 
ignored if 103 is .It. 0., which means a planform 
shape is to be read in. 

The trailing edge sweep in degrees 

Root dimension or chord length along the symmetry axis 
<0. the planform shape is read in from 241-320 
>0. this value is used to calculate the geometry 
The i value at the start of the root. 

The z value for the root. This value will be used 
with 109 or will be added to values from 121-160. 

The value of the SPAN. Used with even spacing. 

If < 0. the chordwise panel spacing is even 

If = 0. the chordwise panel spacing is cosine 

If > 0. the chordwise panel spacing is half cosine 

The x/c for the flap hinge 

The dihedral angle in degrees, (used with 105) 

The spanwise control point location. If 0., the 
centroid is used unless location 201 is nonzero. 

= 0. Symmetry about y=0. is assumed 

* 0. No symmetry 

The number of chordwise locations of camber input. 

The number of spanwise locations of camber input, 
if MD < 3 only the first value will be used, 
i.e. all span stations will have the same camber. 


Planform Shape 

z values at the y coordinates 

y coordinates, (if ( (Y(2)-Y(1))**2+(Z(2)-Z(1) )**2) is 
0. The semi span is broken into M equal segments) 
there must be M+l values input 

y coordinates for the control points. Nonzero values 
will be used to override values based on location 110. 


The leading edge coordinates at Y(J) 

These values are considered only if 103 is < 0 
XLE(l) Corresponds to the coordinate on the axis 
XLE(m+l) corresponds to the coordinate at the tip 
any values which are exactly 0. will be changed 
to make the edge straight between the two 
surrounding nonzero values. 

The trailing edge coordinates. Same format as 241-280 
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321 ZD(I) The x/c values where the camber is input. See 112 

341 YD(J) The y values where the camber is input. See 113 

361 ZD(J) The z values where the camber is input. 

381 TWIST(J) The twist in degrees at the above values of (y,z) 

401 The values of dz/dx for the camber. 

401 to 400+ND are for YD(1) at XD(1) to ZD(ND) 

401+ND to 400+ND*2 are for YD(2) at ZD(1) to ZD(ND) 
These values are interpolated to obtain the boundary 
conditions at the control points. 


Body Geometry 


702 NY 

703 NZ 

704 NZ 


705 NC 


The number of panels around the body (half). 

The number of panels in the x direction on the body. 
The number of singularity segments on the slender 
body. If < 0 the singularity segments will be 
shifted along the Mach lines (supersonic only) 

The number of control points on- the slender body. 


706 IJ.K 


I = 

the order 

of the 

J = 

the order 

of the 

If 

I = 0 

I = 1 

If 

I > 2 

I = 2 

If 

J = 0 

J = 2 

If 

J > 3 

J = 3 


singularities (1,2) 


is used. 


is used, 
is used. 


There is always a 1st order line source at the nose. 
There is always a 2nd order line doublet at the nose. 


IJ < 0 exact conical solution desired at nose. 

IJ > 30 the nose 1st order line source strength = 0. 

the nose 2nd order line doublet strength =0. 

K > 0 a second order axial solution is performed. 

K * 1 no r*phir**3 term included in order 2 solution. 

K = 2 r*phir**3 term included in order 2 solution. 


707 ECC 

708 = 0. 
> 0 . 


The eccentricity of the body cross sections. 

Area = pi * r(I)**2 for body cross sections. 

Velocity boundary conditions on body (beta2x = 1.0) 
Mass flux boundary conditions on body. (beta2x = beta2) 
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709 M.N > 0 Print source and doublet solutions. 

M > 1 Print shifted singularity points and loads. 

> 2 Print source and doublet coefficients. 

= 1 Print source influence matrix 

N = 2 Print doublet influence matrix 

> 2 Print both influence matricies 

710 10 unit number for placing axial loads in a plot dataset 

(use unit 12) 

711 WBI(l) - 0. the body and 1st surface intersection is computed. 

712 WBI (2) - 0. the body and 2nd surface intersection is computed. 

713 WBI (3) = 0. the body and 3rd surface intersection is computed. 

714 WBI (4) = 0. the body and 4th surface intersection is computed. 


715 DELTA Used to calculate the maximum allowable slope of 

bodies. For the axial singularity calculation a 
conical nose extension is constructed tangent to the 
body such that: 

dr/dx = 1.0 / beta - DELTA if delta > 0. 

This operation is performed only for Mach > 1.0 . 
on the region of the actual body where, 

dr/dx > 1.0 / beta - DELTA 

Tangent cone formulas are used to calculate Cp's. 


The type of pressure coefficient calculation, and the 
angle of attach, for 721-726, are determined from 
locations 51-58 and 61-68. 


721 THETA - 1 the 1st theta for Cp computation on body. 

722 THETA - 2 the 2nd theta for Cp computation on body. 

723 THETA - 3 the 3rd theta for Cp computation on body. 

724 THETA - 4 the 4th theta for Cp computation on body. 

725 THETA - 5 the 5th theta for Cp computation on body. 

726 THETA - 6 the 6th theta for Cp computation on body. 
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741 VSHELL(ISF) 

= 1. the surface is a shell composed of vortex panels 
only the upper surface is considered. 

= 0. the surface is a normal lifting surface. 


751 

XINLET(I) 

801 

XG(l ) 

851 

XX(I) 

901 

RR(I) 

1001 

X(I) 

1101 

XC(I) 

1201 

CAM(IJ) 


— ..... ~ as uaauuicu lv nave all 111J.C1 

the mass flow is (l.-XINLET(I)). 

The x coordinates of the body geometry sections. 

The x coordinates of the body panel cross sections, 

if XX(I) = 0., and DATA(703) > 0., XX(I) * XG(I) 
is assumed. 

The r values of the body geometry sections. 

The x coordinates of the slender body singularities. 

The x coordinates of the slender body control points. 

The camber for each panel. Input as dz/dx at the 
control point (each surface input separately), 
unless the input is being read through an APAS 
dataset. For an APAS input all camber values 
are input at one time and the value of DATA( 14) 
must be considered. 


i.e. If the data is not being input through an APAS 

dataset, DATA( 1200+1) is the initial camber dz/dx 
for the I'th panel of the surface being input. 

1801 TWIST(J) The twist for each span station. This is for APAS 
dataset inputs only, and the twist values for all 
span stations are input at one time. See DATA(14) 
for additional information on twist input. 
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1996 FLOW ( 1 ) Mass flow coefficient for inlet # 1 

1997 FLOW ( 2 ) Mass flow coefficient for inlet | 2 

1998 FL0W(3) Mass flow coefficient for inlet # 3 

1999 FLOW ( 4 ) Mass flow coefficient for inlet # 4 

2000 FLOW(5) Mass flow coefficient for inlet # 5 

FLOW(I) * Average value of (l.-beta2x*u) for all field point 
of inlet I (see location 708). 


maximum of 200 field points 

2001 x, y, z, inlet # for field point #1 

2006 x, y, z, inlet f for field point #2 

2011 x f j, z, inlet # for field point #3 

2016 x, y, z, inlet # for field point #4 

2021 x, y f z, inlet # for field point |5 

2026 x, j, z, inlet f for field point #6 

inlet # = 0 means the point is not an inlet point. 
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SUBROUTINE DECRD1 


Subroutine DECRD1 is used to read input data from a fixed block dataset. 
The input data which is read is placed in floating point format into locations 
of the array "DATA" which appears in the argument list. The subroutine will 
also read a single title card of up to 72 characters and place the literal 
data in the array "TITLE" which also appears in the argument list. Subroutine 
DECRD1 allows input in two different formats, called decrd format and free 
format. A description of these input formats will follow. 

The first time DECRD1 is entered, the first record is assumed to be in 
decrd format with LRECL — 72. LRECL is the length of each record which is 
read (maximum LRECL = 132). 

If the characters "C","D", or "F" appear in column 1, the card will not 
be read for any data (except for LRECL). 

A ’C' indicates a comment card. 

A 'D' indicates the following cards are to be read using decrd format. 

An 'F' indicates the following cards are to be read using free format. 

The value of LRECL appears in the card with an F in column 1 in the form 

LRECL(KLM), where KLM is a three digit integer. It will remain constant each 

time DECRD1 is entered, unless it is reset. 

Unless the first card containing data in each entry to DECRD 1 has a blank 
or a minus sign appearing in column 1, with blanks in columns 2-4, the first 
card is assumed to be a title card containing 72 characters of literal data. 
However, cards with the characters "C", "D", or "F" in column 1 do not count 
as data cards and may appear before the title card. When using free format, 
the first data card of an entry to DECRD 1 should be checked carefully to see 
that the first four columns are of the correct form, to avoid confusion 
between a card containing title data and a card with numerical data. 


DECRD FORMAT 

In Decrd format each data card has the format I12,5E12.5. The first 
number on each card is an integer giving a location in the input array; the 
following numbers on the card specify the values to be input into that and the 
four immediately succeeding locations. These numbers can be input either in F 
format, which must include a decimal point or in E12.5 format. Locations left 
blank will remain unchanged. The last card to be read for each call to DECRD1 
is signaled by having a negative location number. 
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FREE FORMAT 


Data may be written in I, F, E, or D format separated by blanks, commas, 
or semicolons. The data will be placed in consecutive locations of the array 
"DATA” in the program. Data locations may be skipped by writing an 1 or using 
consecutive commas with nothing other than blanks between. When data 
locations are skipped, the previous values remain unchanged. 

If the first piece of data on a card is a positive or negative integer, 
which is not a multiplying factor to be described below, it will be assumed to 
be the location number of the next piece of data on the card. If this integer 
is negative this card will be the last card read until subroutine DECRD1 is 
called again. This means the last card must have a negative location number 
as the first piece of data on the card. If the piece of data is not an 
integer, or if it is a multiplying factor, the data will be placed in 
consecutive locations continuing from the previous card. 

A semicolon can be used to designate the end of one card and the start of 
a new card. This means that the next piece of data will be treated as if it 
were the first piece of data on a new card. For this new card everything in 
the previous paragraph applies. 

An integer followed by *■* will cause the subsequent data to be placed in 
consecutive locations beginning with the integer. Blanks may be placed before 
or after the '=’. 

A * preceeded by an integer is used to designate a multiplying factor, 
and will cause subsequent data to be repeated the integer number of times. 

E.g. 5 * 4.4 will result in the value 4.4 being placed in 5 consecutive data 
locations starting with the appropriate data location. Writing 10*X will 
result in 10 data locations being skipped. Spaces before and after a * are 
ignored. A multiplier must be an integer (I format), and it must be > 1. 

Everything within ( ) will be ignored, unless it is preceeded by a 
multiplying factor which is greater than one. E.g. 5*( 3.2 4.6 x 5 ) will 
cause the sequence of data ( 3. 24. 6x5 ) to be repeated 5 times. Nested 
parentheses are not allowed. 
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LIST OF VARIABLES 
GENERAL VARIABLES 


CHORD(J) The chord value at the centroid of span station J 
IJS(IS) The value of IJ where section IS begins 

IJO(J) The number of the first panel at span station J 

ISO(ISF) The section number where surface ISF begins. 

JS(IS) The span station where section IS begins 

MS(IS) The number of span stations in section IS 

NB The number of basic solutions. 

NCHORD(J) The number of panels at span station J 

NS The total number of sections 

NSF The total number of surfaces 

NSPAN(ISF) The number of span stations on surface ISF 

NSL The number of vortex shell sections 

NST The total number of span stations on all lifting 

surfaces. 

NSTS The number of source span stations 

NTB The number of body panels. 

NTL The number of vortex shell panels. 

NTP The total number of panels. 

NTS The number of source parameters. 

NTSL The total number of vortex shell span stations 

NTV The number of vortex panels. 


LIFTING SURFACE PANELS 


CAM(IJ) The normal velocity at the control point due to 

COSZ(IJ) The cosine of the normal of panel IJ 

CX(I r IS) x/c values for the control points of section IS 
DS(J) Width of span station J 

DWDX(IJ) The local dw/dx at thickness control point IJ 
DX(I,IS) Delta (x/c) for the panels of section IS 

ETA(J) Fraction of span running distance for (YCC(J),ZCC(J)) 

camber. 

PA(IJ) The area of panel IJ 

SINZ(IJ) The sine of the normal of panel IJ 

TWIST(J) Twist of span station J 

WTHK(IJ) The local dz/dx at thickness control point IJ 

X(KL,IC) x values of the corners of vortex panel KL (4) 

XC(IJ,l) x value of the centroid for vortex panel IJ 

XC(IJ,2) x value of the control point for vortex panel IJ 

XLE(1,J) x value at the leading left edge of span station J. 

XLE(2,J) x value at the leading right edge of span station J. 

XM(I,IS) x/c values for the midpoints of panels of section IS 

XTE(1,J) x value at the trailing left edge of span station J. 
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XTE(2,J) z value at the trailing right edge of span station J. 

XX(I,IS) x/c values for the panels of section IS 

Y(RL,IC) y values of the corners of vortex panel KL (2) 

YC(IJ) y value of the control point for vortex panel IJ 

YCC(J) y value at the control point of span station J 

YO ( 1 , J) y value at the left edge of span station J. 

Y0(2 f J) y value at the right edge of span station J. 

Z(KL,IC) z values of the corners of vortex panel KL (2) 

ZC(IJ) z value of the control point for vortex panel IJ 

ZCC(J) z value at the control point of span station J 

ZTHK(IJ) The local z/c at thickness control point IJ 
Z0(1,J) z value at the left edge of span station J. 

Z0(2,J) z value at the right edge of span station J. 


BODY PANELS 

B(I) = SQRT(DY(I)**2 + BETA2*DX(l )**2) I = 1,4 
BT(IC,IJ) tans**2 + beta2 for each side of body panel IJ 
= 1.- Mach**2 mass flux boundary conditions. 

i.e. (1. + beta2x *u)*XN+v*YN+(w+ alpha ) * ZN = 0. 

BETA2X ■ 1. velocity boundary conditions. 

IBB(K) The body station number of the first station in body 
section K. 

IJB(K) The panel number of the first body panel in section K. 
NX(K) The number of x body stations for body section K 

NY(K) The number of panels at each body station (half) for 

body section K. 

XB(IC,KL) x values of the corners of body panel KL (4) 

XINLET(IJ) mass flow coefficient of body panel IJ; 0. = impermeable 
XN(1,IJ) » x - component of body panel IJ normal 
XN(2,IJ) * y - component of body panel IJ normal 
XN(3,IJ) - z - component of body panel IJ normal 

XN(4) = XN( 1 ) / XN(8) 

XN(5) = XN(2) / SQRT(XN(2)**2 + XN(3)**2) 

XN(6) = XN(3) / SQRT(XN(2)**2 + XN(3)**2) 

XN(7) = SQRT(XN{2)**2 + XN(3)**2) / XN(8) 

XN(8) = SQRT ( beta2*XN ( 1 ) **2 + XN(2)**2 + XN(3)**2) 

XN(9) = SQRT( XN(1)**2 + XN(2)**2 + XN(3)**2) = 2. * area 

XP(IC,KL) x values of the corners of body panel KL (4) planar 

XBO(I) The x value of the center of body station I. 

X0(1,IJ) x value of the centroid of body panel IJ 
X0(2,IJ) y value of the centroid of body panel IJ 
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X0(3,IJ) z value of the centroid of body panel IJ 
YB(IC,KL) x values of the corners of body panel KL (4) 

YP(IC,KL) x values of the corners of body panel KL (4) planar 

ZB(IC,KL) x values of the corners of body panel KL (4) 

ZP(IC,KL) x values of the corners of body panel KL (4) planar 


COMPUTED VARIABLES 

Cp(IJ,K) The delta-Cp across panel IJ from basic solution K. 

UK(IJ,R) The control point upper surface z velocity. 

VK(IJ,K) The control point upper surface binormal velocity. 
WK(IJ,K) The control point upper surface normal velocity. 

PR(IJ,K) The control point upper surface velocity potential. 

US(IJ) The thickness induced upper surface z velocity. 

VS(IJ) The thickness induced upper surface binormal velocity. 

WS(IJ) The thickness induced upper surface normal velocity. 

PS(IJ) The thickness induced upper surface velocity potential. 

The following variables are required for 2nd order solutions only. 

WSX(IJ) d/dz of WS(IJ) source normal velocity. - 

WCX(IJ) d/dz of WK(IJ,2) camber normal velocity. 

KKX The number of 2nd order boundary condition solutions. 

* 6 if there is no body. 

=9 if there is a body. 

UBE(IJ,K) The even symmetry u velocities due to 2nd order b.c. 
UBO(IJ,K) The odd symmetry u velocities due to 2nd order b.c. 
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SUBPROGRAMS 


MAIN 

ABMAIN 

ACMAIN 

AER02 

AER02B 

AER02P 

AER02V 

AMINMX 

APAS 

APASB 

AXXABA 

AXXAIJ 

AZZBDT 

AZZBXR 

AXXCSD 

AXXLIN 

AXXLOD 

AXXNCP 

AXXNFO 

AXXNF1 

AXXNOS 


AXXSLX 

AXXTCN 

AXXUVW 

AXXVEL 

AXXZIS 

BDTAIJ 

BDYSCE 


FUNCTION 

Sets the main array size for the program. 

Reads input data, computes geometry, and sets array 
sizes. 

Controls the flow of the program. 

Controls the program flow for the computation of 1st 
or 2nd order pressures and loads. 

Computes the u velocities due to 2nd order boundary 
Conditions. 

Computes the 1st or 2nd order pressures from the 
induced velocities. 

Computes the odd and even symmetry velocities induced 
by the basic solutions. 

Finds the largest, smallest or largest absolute value 
of the elements of an array. 

Reads geometry from APAS output. 

Modifies APAS body geometry coordinates to form desired 
panel coordinates. 

Checks if an axis singularity solution is desired. 

Computes axis singularity source and doublet 
influence coefficients. 

Computes body geometry for axis singularity solution. 

Prints axisymmetric body geometry and coefficients of 
source and doublet coefficients. 

Computes coefficients of semi-infinite source and 
doublet line singularities. 

Computes x and r for axisymmetric body geometry. 

Computes Cp's and forces on isolated body from axis 
singularity solution. 

Integrates forces on body nose using Cp's from a 
tangent cone solution. 

Function for body nose curve fit. 

Function for body nose curve fit. 

Finds the point where a conical nose extension is 
added to the body in order to avoid body slopes in 
excess of the Mach angle. 

Controls program flow for axis singularity solution. 

Computes tangent cone Cp's on the body. 

Computes velocities from axis singularities. 

Computes velocities at panel control points from 
axis singularities. 

Computes influence coefficients for source and doublet 
semi-infinite line singularities. 

Computes influence coefficients for body panels. 

Calculates u,v,w,phi influence coefficients for a 
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B0D7DS 

B0D7G 

B0D7GM 

BODYRD 

BODYW 

CLCM 

DECRD1 

DISPLY 

DMXMVE 

DSSPLY 

ENDREC 

FIELD 


FXDX3 

GEOM 

GOMTRY 

HSHLDR 

INLET 

INTRP3 

INTRP4 

INTRPX 

INTRPY 

INTSCT 

LINEAR 

MATRIX 

MATRXF 

MATRXT 

MTXADD 

MTXMLT 

MTXMVE 


unit strength source panel. 

Prints arrays of characteristics at body panel control 
points. 

Computes body panel geometric characteristics. 

Computes and prints body geometry from input data. 

Reads body panel geometry from APAS dataset. 

Finds the intersection of the body and the aerodynmaic 
surfaces. 

Computes the lift drag and moment characteristics of 
aerodynamic surfaces. 

Reads the input data. 

Prints arrays of characteristics at the panel control 
points of the aerodynmaic surfaces. 

Moves the elements of a double precision array. 

Prints arrays of characteristics at the panel corner 
points of the aerodynamic surfaces. 

Used with errset to check for APAS type influence 
matrix. 

Computes first order field point properties using 
previously calculated influence coefficients and 
the first order solution. 

Integrates data using a third order curve fit through 
the nearest four points. 

Computes the geometric characteristics of the 
aerodynamic surfaces and panels. 

Computes the panel geometry for the aerodynamic 

- surfaces. 

Solves sets of linear simultaneous equations in a 
least square sense. 

Performs solution while satisfying boundary conditions 
at a set of inlet points. 

Interpolates or differentiates data using a third order 
curve fit through the nearest four points. 

Interpolates or differentiates data using a third or 
fourth order curve fit through the nearest four points. 

Interpolates or differentiates properties chordwise on an 
aerodynamic surface using subroutine INTRP4. 

Interpolates or differentiates properties spanwise on an 
aerodynamic surface using subroutine INTRP4. 

Finds the intersection of two lines. 

Fills in aerodynamic surface geometry not input, and 
Calculates geometry. 

Computes aerodynamic influence coefficients for 
source and vortex panels on aerodynamic surfaces. 

Displays arrays of data. 

Displays arrays of data. 

Adds multiples of two arrays. 

Multiplies two arrays. 

Moves the elements of an array. 
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NOTZRO 

ORTHO 


PHIS11 

PHIXY 

PLOT 

POLAR 

QUADPL 

SQFIT 

THICK 

VORTEX 

VTXAIJ 

VTXDRG 

VTXDR2 

VTXLFT 

HBINF 

WBSOL 

WBUVW 

WBVEL 

WINTRP 


file 

5 

6 
8 
9 

10 

11 

12 

13 


Checks to see if an array has any nonzero elements. 

Solves sets of linear simultaneous equations using the 
method of succesive orthogonalizations. A 
quasi-inverse matrix may be computed or, if 
previously computed, used to perform the solution. 

Calculates u,v,w,phi influence coefficients for aero 
surface source panel edges on control points. 

Calculates v (binormal) velocities on aerodynamic 
surfaces from phi (for 2nd order solution). 

Writes geometry and aerodynamic data on a disk unit 
for computer graphics. 

Calculates first order drag polar and aerodynamic 
cofficients. 

Computes body source panel geometric characteristics. 

Calculates leading edge suction using a least square 
Curve fit of Cp to cot(x/c). 

Computes wing thickness z/c, dz/dx, d2z/dx2 from input. 

Calculates u,v,w,phi influence coefficients for aero 
surface vortex panel edges on control points. 

Controls program flow for aerodynamic influence 
Coefficient calculation. 

Calculates vortex drag in the Trefftz plane. 

Calculates coefficients for VTXDRG. 

Calculates vortex lift from leading edge suction. 

Prints source or vortex influence matricies for 
aerodynamic surfaces. 

Controls the matrix solution using subroutine ortho. 

Calculates boundary conditions for solution. 

Calculates u,v,v,phi at panel control points. 

Interpolates velocities to various points on panels. 


use 

input data 
Print output 
scratch file 
scratch file 
scratch file 

storage of influence matricies and quasi-inverse 

plot output file 

APAS geometry input file 
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FLOW DIAGRAM 


********** 

* AAMAIN *' 

********** 

I 

V 

********** 

* ABMAIN *<« 
********** 


********** 

<==>* decrdi * 

********** 

F 

********** 

<==>* APAS *<> 

********** 

I 

********** 

<==>* LINEAR * 

********** 

r 

********** 

<==>* GOMTRY * 

********** 

F 

I ********** 

<==>* THICK * 

********** 


********** 

===>* APASB * 
********** 
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V 

********** 

* ACMAIN *<==> 

********** 


********** 

<==>* GEOM * 

********** 


********** 

<==>* bodyg *< 

********** 


********** 

==>)<==>* BODYRD * 

********** 


********** 

<==>* AXXLIN * 

********** 

f 

********** 

<==>* QUADPL * 

********** 

r 

********** 

<==>* BODYGM * 

********** 

f 

********** 

<==>* BODYW * 

********** 

********** 

<==>* AXXABA *<==> see Slender Body Flow diagram 

********** 


********** 

<==>* VTXAIJ *<= 
********** 


********** ********** 

=>* MATRIX *<==>!<==>* VORTEX * 


********** 


| ********** 

********** 

<==>* PHI SI 1 * 

********** 
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********** ********** 

<==>* BDTAIJ *<=======>* BDYSCE * 

********** ********** 

V 

********** 

<==>* WBUVW * 

********** 

V 

********** ********** 

<==>* WBSOL *<sa«K»>* ORTHO * 

********** ********** 

V 

********** 

<=■>* WBVEL * 

********** 

V 

********** 

<==>* WBINF * 

********** 

V 

********** ********** 

<==>* VTXDRG *<=======>* VTXDR2 * 

********** ********** 

V 
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********** 

<==>* AER02 *< 

********** 

V 


V 


V 


V 


V 


V 

********** 

<==>* POLAR * 

********** 


********** 

> <==>* PHIXY * 

********** 

V 

********** 

<==>* AER02V * 

********** 

V 

********** 

>* WBSOL * 

********** 

********** 

* ORTHO * 

********** 

V 

********** ********** 

<==>* CLCM *<=>[<=>* SQFIT * 

********** ********** 

V V 

I ********** I ********** 

<==>* PLOT * |<=>* VTXLFT * 

********** ********** 


********** 

<==>* AER02B *< 

********** 

V 

********** 

<«=>* AER02P * 

********** 
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SLENDER BODY CALCULATION 


********** 

* AXXABA * 

********** 


********** 

* AXXBDY *<==> 

********** 


********** 

<==>* AXXLIN * 

********** 

r 

********** 

<==>* AXXNOS *<=== 

********** 

r 

********** 

<==>* AXXSLX *<==> 

********** 


********** 

<==>* AXXVEL *<=== 

********** 


********** 

===>* AXXNCP *<= 

********** 

' ********** 

<==>* AXXCSD * 

********** 

r 

********** 

<==>* AXXBXR * 

********** 

r 

********** 

<«>* AXXAIJ *<= 

********** 

r 

********** 

<==>* HSHLDR * 

********** 

r 

********** 

<==>* AXXUVW * 

********** 

f 

I ********** 

<==>* AXXLOD * 

********** 


********** 

====>* AXXAIJ *<=> 

********** 


********** 

>* AXXTCN * 

********** 


********** 

>* AXXXIS * 

********** 


********** 

<=>* AXXXIS * 

********** 

********** 

<=>* AXXCOR * 

********** 
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TEST CASE 


Results for the seventy degree delta wing-body arrangement of figure 1 
are presented in this section. 

The aerodynamic paneling representation used in the analysis is presented 
on figure 2 and is typical of the program geometrical graphics output. 

Comparison of first and second order results with experiments for a Mach 
number of 6 and an angle of attack of 8 degrees are presented, on figure 3. 
Improved wing surface pressure coefficient predictions are systematically 
obtained for the second order analysis with the exception of the root section 
on the compression side and are in reasonably good agreement with 
measurements. Additional results are presented in references 2 and 3. 

Test case input is presented on pages 41 through 43, and detailed program 
printed output is presented on pages 44 through 98. Typical aerodynamic data 
graphics output is presented on pages 99 through 109. 
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Figure 2. Wing-Body Finite Element Model 
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Figure 3. Continued 

















Figure 3. Continued 
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TEST CASE INPUT DATA 


WBODT . DATA { M600 ) 1ST-AX (10 X 11) (8 X 11) VEL B.C. NEG 37X37 
C CALCULATE NEW AERO MATRIX AND INVERSE AND SAVE (=2.0) 

7 2. 

C IF < 0. CALCULATES A SECOND ORDER SOLUTION. 

8 - 1 . 

C = 2.0 FOR SPANWISE AND CHORDWISE LINEARLY VARING SOURCE PANELS. 

10 2 . 

C GEOMETRY CALCULATION ONLY 

12 0 . 

C CONTROL POINT AT PANEL CENTROID (VORTEX PANELS) 

13 1. 

C NO PRINTING OF SOURCE (THICKNESS) DZDX AND Z/C MATRICIES 
15 0. 

=3.0 PRINTS ALL PERTURBATION VELOCITIES DUE TO THICKNESS 

18 3.5 

PRINTS ALMOST EVERYTHING (EACH LOWER NUMBER PRINTS LESS) 

19 0. 

PLOT DATASET (EXTENDED) NO PLOT DATASET (= 0.) 

20 3. 

EXTENDS FIRST SURFACE TO PANEL CENTERLINE WITH NEGATIVE SWEEP 
26 1.2 
MACH NUMBER 
35 6.0 

PRINTS ODD AND EVEN SYMMETRY VELOCITIES 

45 334. 

46 334. 

47 334. 

PRINTS CAMBER AND THICKNESS NORMAL VELOCITIES 

48 0.0 

ANGLE OF ATTACK ( > 90. FOR ADD LOAD ) 

51 99.0 99.0 8.0 8.0 

1ST ORDER AXIAL SOLUTION ( > 0. ) 

ISENTROPIC =4 (ON BODY) 

1ST ORDER = 1 (ON LIFTING SURFACES) 

NO CAMBER = 0 

NO THICKNESS = 0 

61 41.00 42.00 41.00 42.00 
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noon 


ORIGINAL PAGE IS 
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C NUMBER OF PANELS AROUND THE BODY (HALF) 

702 8. 

AXIAL SINGULARITIES. 

< 0 MEANS EXACT CONICAL SOLUTION AT NOSE. 

2 = AXIAL VARIATION ORDER OF SOURCE SINGULARITIES 

3 = AXIAL VARIATION ORDER OF DOUBLET SINGULARITIES 

706-23.0 
C VEL B.C. = 0. 

708 0. 

C PRINT SOURCE AND DOUBLET SOLUTIONS. 

709 0. 

710 0.0 

C IF > 0. NO BODY AND SURFACE INTERSECTION IS COMPUTED 

711 1. 1. 1. 1. 1. 

C MAXIMUM ALLOWABLE SLOPE ON BODY * 1. / BETA - DATA(715) 

715 0.0029 

C PRESSURE COEFFICIENT CALCULATION ON BODY (AXIAL SOLUTION) 
720 4.0 

C ANGLE FOR PRESSURE COEFFICIENT CALCULATION. 



721 0. 

22.5 

45. 

67.5 

90. 

c 

X VALUES OF THE 1 

BODY GEOMETRY 

SECTIONS. 




801 0.0 

1.00 

1.09 

1.10 

1.20 


806 1.30 

1.40 

1.50 

2.00 

3.00 


811 4.00 

5.00 

6.00 

6.70 

6.75 


816 6.80 

6.810 

6.820 

11.10 

11.150 


821 11.203 

11.30 

11.50 

12.00 

12.50 


826 13.00 

13.50 

13.983 

14.10 

14.20 


831 14.30 

14.40 

14.50 

14.60 

14.70 


836 14.80 

16.9360 




c 

X COORDINATES OF 

THE BODY PANEL CROSS SECTIONS (R INTERPOLATED) 


851 6.1839 






862 14.0775 





c 

R VALUES OF THE i 

BODY GEOMETRY 

SECTIONS. 




901 0. 

0.250501 

0.265452 

0.267113 

0.283416 


906 0.299132 

0.314302 

0.328960 

0.395480 

0.501051 


911 0.577263 

0.628500 

0.656942 

0.663870 

0.663961 


916 0.663999 

0.664 

0.664 

0.664 

0.664 


921 0.664 

0.663791 

0.662036 

0.649816 

0.626208 


926 0.590771 

0.542754 

0.483287 

0.466763 

0.451935 


931 0.436431 

0.420224 

0.403653 

0.387083 
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WING 

C CHORDWISE SPANWISE # OP PANELS 

2 10 . 11 . 

C T/C A THICKNESS MATRIX IS CALCULATED (UNLESS * 0.) 

89 0.05 -1.0 

C ROOT ( < 0. INDICATES LEADING AND TRAILING EDGE INPUT) 

103-9.00 

C <0. INDICATE PANEL SPACING IS EVEN. 

107-1.0 

C T CORDI NATES OF PANEL ENDS 
161 0.6640 0.9323 

172 3.2215 

C X CORDI NATES OP LEADING EDGE 
241 6.8072 

252 13.8333 

C X CORDINATES OP TRAILING EDGE 
281 13.983 
292 13.983 
292 13.983 

C >1.0 LAST SURFACE HAS BEEN READ. 

1 1.0 

C - -1.0 NO MORE CASES FOLLOW. 

1 - 1 . 
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ORIGINAL 
OF POOR QUAU 


DATA SET 4 
MACH *6.20 
3 * 3. a 

AL°H4 a 8.00 DEG 
NS?C-«»*M*aMa« 


P-OT PILE created BY OPT 


ENTER! j to plot 
0 TO BYPASS this DATA SET 


ENTER • TO PROCEED TO THE NEXT CASE. 


25 2-9J Ef TMICXNESS distribution. 

S 2-OT THE CAMBER DISTRIBUTION. 

TO PLOT Z/C'S FROM TWIST V CAMBER. 

S PLOT OpScT^ 0 ^ TWIST C4MBER 4 FLAPS * 
TO PLOT CP UPPER. 

TO PLOT CP BOTH. 

2 2-2J £ upper cp lower 
TO PLOT SPANWISE CHARACTERISTICS. 


> 8 TO EXIT CLIST. 


? 

8 
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WBOD v . DATA( K800 ) 1ST-AX ( 10 X 1 1 ) (8 X 1 1 ) VEL B.C. NE6 37X37 

ALPHA ■ 8.000 CP-UP AND CP-lO STATION * I 

I 0 2. 252 



100 









WBOD v . DATA( M600 ) 1ST-AX (10 X !1) (8 X II) VEL B.C 

ALPHA ■ 0.000 CP-UP and CP-lO S’A'ICN • 3 

3 0 0.238 


ORIGINAL RASE ig 
OR POOR QUALITY 
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0.140 





















WBODY . DA"A( K603 ) 1ST-4X (13 X 11) (8 X 1!) VEL B.C 
aipha B.oaa cp-cp an: cp-lo s t a*icn • 7 








ORIGINAL PAGH is. 

OF POOR QUALITY 


I 

I 

I 








• TO PROCEED TO Tt€ NEXT CASE 





WBODY. DATA (1^600) 1ST-AX (10 X 11M8 X 1:) VEL B.C 

A..PHA = 8.000 NORMAL FORCE S-R C ACE * 1 









WBODY.DAT A (M600) 1ST-AX (10 X 11) (8 X 1!) VEL B.C 

<VP.iA * 8.080 lOCAu CL SURFACE « 1 



0.203 

















MBO0) 1ST-AX (10 X 11) (8 X 1!) VEL B.C. NEG 37X37 

SPAWISE WAS SURFACE • 1 
















OPTIMIZATION PROGRAM OPT 


DISCUSSION 


The computer program OPT obtains the optimum camber, twist, and/or flap 
deflections for a given configuration by obtaining the minimum zero suction 
drag subject to various constraints. On the lifting surface panels, the 
component of force in the drag direction is equal to the local angle of 
incidence times the force in the direction normal to the panel. The zero 
suction drag is thus defined as the sum: 


N 

Cd * SREF = sum Cp(i) * Alpha(i) * Pa(i) 
i=l 

The index i runs over all N finite elements (panels) into which the 
aerodynamic surfaces are divided, and 

Cp(i) Is the delta Cp across element i 
Alpha(i) Is the local angle of attack of element i 
Pa(i) Is the area of element i 

All pressures and velocities in program OPT are computed using 
subroutines from program WBODY which, for first order analysis, uses all of 
the assumptions of linearized, thin wing theory. Optimizations may be 
performed at all Mach numbers on configurations having up to ten surfaces, and 
one body. The program does not compute configuration geometry or aerodynamic 
influence matricies. This information is read from data files which are 
created by other programs such as VBODY or APAS. 

The final optimized result may be constrained to maintain a total lift 
coefficient and center of pressure, and various combinations of lift 
coefficients on any specified set of surfaces, spanwise section lift 
coefficients, spanwise variation of sectional center of pressure (x/c-C.P.), 
spanwise twist variation, twist and camber over specified span stations, and 
flap deflection angle. 


.... Por a first order optimization, with or without a paneled body, two 
different methods are used. The first uses constraint functions for chordwise 
distributions at any or all span stations. Span stations having no 
chordwise pressure constraint functions may be thought of as actually having a 
number of constraint functions equal to the number of chordwise panels. Each 
function is then a delta function which is equal to 1.0 on the panel equal to 
the constraint function number, and equal to 0.0 elsewhere. 

The second method uses chordwise polynomial curves, and flap deflection 
slopes over selected panels as camber constraint functions. The solution in 
each case consists of the coefficients of these constraint functions at each 
span station. If a paneled body is present, since its' geometrical shape is 
fixed, the only degree of freedom allowed for it is the angle of attack, which 
may be constrained to a fixed value if desired. 

Second order optimizations may be performed, using camber constraint 
functions on planar configurations without a body. The influence of thickness 
must be considered, otherwise the result will be the same as a first order 
result. 

Cases using chordwise pressure constraint functions, or using no 
constraint functions, will be referred to as a Cp optimization. This method 
of solution first yields panel strengths (delta-Cp's on vortex panels and 
source strengths on body source panels) which result in minimum drag, subject 
to any constraints, and then the camber and twist is determined from the panel 
strengths and the aerodynamic influence matrix. 

When using camber constraint' functions, the optimization will be called a 
twist optimization, which is a shortened form of saying twist and camber 
optimization. This name originates because the solution yields the twist and 
camber shapes which result in minimum drag, subject to any constraints. The 
panel strengths (delta-Cp's on vortex panels and source strengths on body 
source panels) may be obtained by solving the set of simultaneous equations 
using the aerodynamic influence matrix. The simplest case using this method 
would be a solution which considered only the lowest order camber functions, 
consisting of uncambered constraint functions (first order polynomials) at 
each span station. This would mean the angle of attack or twist at each span 
station was the only degree of freedom, and the solution would consist of the 
optimum distribution of twist with whatever initial camber was input. 

The program uses geometry data, aerodynamic influence matricies, and the 
quasi-inverse matrix (if present), from either computer program WBODY, APAS or 
TOP. These data are read from datasets which are allocated to the correct 
unit numbers. 

The program may also be used for analysis only using the geometry and 
aerodynamic influence matricies read from the respective data sets. 


Ill 
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INPUT DATA 


Input data, which directs the program operation, is read using subroutine 
DECRD1 , described on page 19, and is stored in the array called "DATA" . All 
locations are initially set equal to 0. Jet flap and vortex lift capability 
are not available in the program at this time and any input which refers to 
these capabilities could cause unpredictable results. 

There are three different types of program operation, a Cp optimization, 
a twist optimization, and an analysis. The following is a brief description 
of the program input necessary for these three different modes. A more 
detailed description of all input locations follows. 


Cp OPTIMIZATION 

This option may be used to find optimum twist and camber subject to 
various constraints. Span stations may be designated where a specified camber 
is to be maintained, or where a specified chordwise Cp distribution is to be 
maintained. The unknowns consist of panel Cp's, at span stations where there 
are no constraint functions, and the coefficients of chordwise Cp functions at 
span stations where there are constraint functions. The camber is obtained by 
multiplying the optimized Cp distribution by the aerodynamic influence matrix. 

The following data must be input for a Cp optimization. 


Input variables 
Total CL 

ZCG - the desired center of pressure 
The following input is optional 


location 

4 

32 


Input variables loc 

The number of chordwise Cp constraint functions 
The type of aero matrix and inverse 
Second order anlysis of optimized result 
Roll, yaw, and side force constraints (asymmetric) 
Reference span station for angle of attack 
Input camber 

Mach number for suction calculation 
x/c for n-velocity in drag analysis & optimization 36 
CBAR 
CAVG 
SREF 


c 

o 

•H 

-M 

default 

5 

2 

7 

0. 

8 

none 

11 

none 

21 

1 

29 

none 

35 

geometry 

36 

0.515,0.875 

40 

geometry 

41 

geometry 

42 

geometry 
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SPAN 

Additional angles of attack 
Additional angles of attack (unit solutions) 
Individual surface CL’s 
Fixed camber, or fixed Cp span stations 
Desired CL*c/Cavg at span stations 
Relative delta CL at span stations 
x/c of center of pressure (C.P.) at span stations 251 
Relative delta x/c-C.P. at span stations 
Relative delta CL at span stations 
Span stations with no constraint functions 
x/c's to input camber (method # 3, see loc 29) 

z/c’s to input camber (method # 3, see loc 29) 

Twist constraint reference stations 
Additional input camber (method # 2, see loc 29) 
or fixed chordwise Cp distributions 
Flap deflections (input camber) 

The default "geometry" means the values from the 
dataset which contain the geometry are used. 


43 


geometry 

52 


none 

62 


none 

71 


none 

81 

+ 

none 

151 

+ 

none 

201 

+ 

none 

; 251 

+ 

none 

301 


none 

351 

+ 

none 

401 

+ 

see 5 

551 

+ 

none 

601 

+ 

none 

951 

+ 

none 

2840 

+ 

none 

3340 

+ 

none 


The following input data controls the output which is printed 


Type of Printout 

Constraint equation printout 

Thickness z/c, dz/dx, and initial camber 

Geometry 

Additional printout 
Plot dataset created 


location default 

14 none 

15 none 

16 some 

19 some 

20 none 


TWIST OPTIMIZATION 

This option may be used to find optimum twist camber and flap 
deflections. The constraint functions consist of camber geometry shapes, 
which result in a series of unit solutions, and the optimization consists of 
finding the coefficients of these unit solutions. The Cp's are obtained by 
solving a set of simultaneous equations using the aerodynamic influence 
matrix. 

The final wing camber will be equal to the input wing camber (see data 
location 29) plus the appropriate sum of the camber constraint functions (if 
any) plus any flap defections. 
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The following data must be input for a twist optimization. 


Input variables location 

Total CL 4 

Twist optimization flag >0. 17 

XCG - the desired center of pressure 32 

The following input is optional 

Input variables location 

The number of camber constraint functions 5 

The type of aero matrix and inverse 7 

Second order optimization or analysis 8 

Roll, yaw, and side force constraints (asymmetric) 11 
Trimmed or not trimmed 12 

Reference span station for angle of attack 21 

Type of interpolation for flaps 23 

Input camber 29 

Mach number for suction calculation 35 

x/c for n-velocity in drag calculation and OPT 36 

CBAR 40 

CAVG 41 

SREF 42 

SPAN 43 

Additional angles of attack 52 

Additional angles of attack (unit solutions) 62 

Individual surface CL's 71 

Flap panel locations 81 + 

Constraints on camber functions 151 + 

Constraints on twist 351 + 

Zero input camber at span stations 401 + 

x/c's to input camber (method I 3, see loc 29) 551 + 

z/c's to input camber (method f 3, see loc 29) 601 + 

Twist constraint reference stations 951 + 

Additional input camber (method # 2, see loc 29) 2840 + 

Flap panel ratios and deflections 3340 + 


The default "geometry" means the values from the 
dataset which contain the geometry are used. 


default 

1 

0 . 

none 

none 

not 

1 

none 

none 

geometry 

0.515,0.875 

geometry 

geometry 

geometry 

geometry 

none 

none 

none 

none 

none 

none 

none 

none 

none 

none 

1.0 

none 


The print control input is the same as for a Cp optimization. 
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ANALYSIS 


The program may be used to analyze a given configuration at various 
angles of attack with or without camber. The program uses the specified input 
geometry and aerodynamic influence matrix. The Cp's at specified span 
stations may be set equal to zero in order to find the influence of various 
span stations or surfaces on other span stations or surfaces. 


The following data input is required for an analysis. 

Input variable location 

Analysis flag ■ 1.0 12 

The following input is optional 
Input variables 


location default 


7 

8 
11 
21 
23 
29 

35 

36 

40 

41 

42 

43 
51 


The type of aero matrix and inverse 
Second order analysis 

Roll, yaw, and side force constraints (asymmetric) 
Reference span station for angle of attack 
Type of interpolation for flaps 
Input camber 

Mach number for suction calculation 
x/c for n-velocity in drag calculation 
CBAR 
CAVG 
SREF 
Span 

Angle of attack 
Unit solutions 
Surface CL's 

Span stations to zero Cp or camber 
Additional input twist 
x/c's to input camber (method # 3, 

z/c's to input camber (method i 3, 

Twist constraint reference stations 
Additional input camber (method # 2, 

Flap panel ratios and deflections 

The default ^geometry" means the values from the 
dataset which contain the geometry are used. 

The print control input is the same as for a Cp optimization. 


0 . 

none 

none 

1 

none 

none 

geometry 
0.515 
geometry 
geometry 
geometry 
geometry 
add load 



61 

41.02 


71 

none 


81 + 

none 


351 + 

none 

see 29) 

551 + 

none 

see 29) 

601 + 

none 


951 + 

none 

see 29) 

2840 + 

1.0 


3340 + 

none 
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DETAILED DESCRIPTION FOR DATA INPUT 


Location Data 

1 If this value is nonzero the case is terminated 

4 The total CL of all surfaces. 

5 I.J NFX = I = The number of chordvise constraint functions 

NFY * J * The number of spanvise constraint functions 

Two cases: 

1. DATA (17) ■ 0. Called a Cp optimization. 

Cp optimized - camber results from Cp 

1=0 NFX * 2 is used as a default value. 

I < 0 No constraint functions are used. 

Constraint functions should always be used for 
subsonic Mach numbers. 

Individual span stations may be specified where no 
constraint functions are desired using data 
locations beginning at 401. 

Fixed span stations (see location 81) always have 
no constraint functions. 

NFY is ignored (NFY = 0 is used regardless of input) 

2. DATA(17) > 0. Called a twist optimization. 

Camber optimized - Cp results from camber 

1=0 Twist only is optimized, (same as I * 1) 

I > 1 Twist + (1-1) camber constraint are used. 

If J * NFY > 0 (works only with a single surface) 

The twist (or camber constraint coefficient) values, 
u(k,eta), are constrained spanwise on each surface 
In the form: 

NFY 

U(k,eta) = sum a(k,i+l) * eta**i k = 1,NFX 
i*0 
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7 


I MATRX. JMATRX 


8 


9 


I MATRX < 0 
I MATRX < - 1 
IMATRZ < - 2 


If the aero matrix is from WB0D7. 

If the last aero matrix is an inverse matrix. 
If 1st aero matrix is a thickness matrix. 


JMATRX > 0 This is used only with configurations having 

bodys which are paneled with source panels. 


JMATRX 

Bodymx 

Body Cp 

Un (for body configurations only) 

0 

Yes 

No 

No 

1 

Yes 

No 

Yes 

2 

Yes 

Yes 

Yes 

3-4 

Half 

No 

No (does constraint functions) 

5 

NO 

No 

No 

6 

No 

Yes 

No 

7-9 

No 

No 

No 

< 0 . 

Second 

order analysis. 


^ Second order optimization is performed. 

< “1. Second order optimization, including thickness drag. 

A second order optimization can be used with DATA( 17) > 0. only 

The value of q (dynamic pressure). This input is used to 
signal that a flexible configuration is being considered. 

This can be used only with a twist optimization 
(DATA(17) >0.). 


10 The number of g’s at the design CL (DATA(4)). This value is 

used to multiply the weights from locations 1801-2301 in 
order to find the deflections due to the inertial loads on 
a flexible configuration. 

11 If 0. The symmetry data from the geometry dataset is ignored. 

=1. Constraint equation included for roll. 

* 2. Constraint equations included for roll and yaw. 

* 3. Constraint equations included for roll yaw and side force. 

>3. Oses symmetry data with no additional constraints. 

These constraints are used only for asymmetric configurations 
XCG, YCG and 2CG obtained from DATA( 32-34) 

If a negative number is input, a solution is obtained for 
flap settings to satisfy a roll constraint for a flexible 
configuration. When applied to symmetric configurations 
an anti-symmetric aero matrix should be used. If q * 0. is 
desired (DATA(9)), set q < 0. 
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12 


If < - 2. 
If = - 2. 
If = - 1. 
If * 0. 
If = 1. 
If > 1. 


Prints geometry only, no optimization or analysis. 
No add load and no drag polar. 

If DATA(17) > 0. , calculates trimmed polar. 

An add load and an optimization is computed. 
Computation using input cambers (analysis only). 
Computation using Cp's from plot dataset (WBODY). 


14 M1X.NUX MIX * the number of constraint equations printed 

NUX = the number of Q(I,J) (drag matrix) rows printed. 
If greater than the number of equations, all are printed. 

If a negative number is input the equations will be printed 
Only if the matrix is singular. 


15 


1. Print source dz/dx and z/c matricies 

2. Print both source and camber arrays 

3. Print camber slope matrix 


16 The geometry printout 

If > 0. or < -1. A wide range of geometric parameters will 
be printed. 


17 If > 0. A twist optimization will be performed, 
unless DATA(12) > 0. 


It may be recalled that a twist optimization is only a 
term used to designate the type of optimization procedure 
involved. The twist optimization procedure may be used to 
find optimum twist camber and flap deflections. The 
following input should be considered when performing a 
twist optimization. 

The final wing camber will be equal to the input wing camber 
(see discussion for data location 29) plus the appropriate 
sum of the camber constraint functions (if any) plus any 
flap defections. 
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18 Printout of velocities due to thickness (source panels). 

* 1. printout of v velocities (normal to panel). 

* 2. printout of u velocities (x direction). 

* 3. printout of v velocities (bi-normal to panel). 

£ 4. printout of z/c for thickness. 

£ 5. printout of dv/dx (2nd order only). 

19 INTMED Various orders of intermediate printout (-2. to 4.) 


20 


-2. Least printout (no drag n-velocities) 

-1. No upper and lover surface cps 

2. More printout 

Prints flexible unit solutions and 
(v-v), (v-b), (b-w), (b-b) 

3. Above + 

4. Above + 2nd order boundary condition solutions 


□PLOT 0. No data is placed in a plot dataset. 

1. Data is placed in an APAS plot dataset (geometry). 

2. Data is placed in a UDP plot dataset. 

> 2. Data is placed in a UDP plot dataset (extended). 


21 The reference span station which will be used to calculate 
the angle of attack. This angle of attack is only used 
for reference purposes in order to calculate twist. If 
A0(K) and A0(J) are the local angles of attack with 
respect to the freest ream, COSZ(J) and COSZ(K) are the local 
dihedral angles, and R is the reference span station, then: 

Alpha * A0(K) / COSZ(R) 

and the value of twist for span station J will be: 

TWIST(J) « A0(J) - Alpha * COSZ(J) 

K < 0 All twist values will be referenced to the freestream, 
i.e. Alpha ■ 0. 

K * 0 R ■ 1 is used. If a body is present then Alpha is 
determined by the angle of attack of the body. 

R > 0 This station is used as the reference span station, 
determined by the angle of attack of the body. 

24 < 0. All normal velocities (slopes), used to compute the zero 

suction drag, are interpolated chordwise from the 
boundary condition control points (usually x/c = 0.875 
of each panel) to x/c » RD (see DATA(36)). A third 
order curve fit is used for this purpose. 
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Otherwise (default), the normal velocities are interpolated 
assuming possible discontinuities where flap panels 
occur (see locations 81-150), and the chordwise 
interpolation of normal velocities is performed while 
the flap deflections are removed. The interpolation is 
performed using information obtained in data location 
36. 


26-28 S.K Surface to be extended to centerline (paneled bodies only) 

If > 0 The surface will only be extended when the aero 
influence matrix is calculated. 

S = surface number (in ascending order). 

R - 0 Extend inboard panels with zero sweep. 

K * 1 Extend inboard panels with same sweep. 

K = 2 Extend inboard panels with negative sweep. 

i.e. if DATA (26) * 2.0 The calculation of the aero 
influence matrix for surface f 2 will be 
performed by extending the inboard panels to 
y ■ 0. with zero sweep. 

29 N.MF (N,M,F integers) Initial camber input specification. 

N < 0 No initial camber is assumed. 

There are three ways of inputing any initial camber and 
the resulting camber slopes from each method will be 
added together. 

1. Reading span station twists and camber normal 

velocities at panel control points from the 
dataset containing the geometry. 

2. Inputing normal velocities at panel control points 

using locations beginning at 2840 or 3340. For 
cases where it is appropriate, values of twist 
may be input using data locations beginning at 
351. 

3. Inputing x/c and z/c using locations beginning at 

3951 and 4001. 
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method # 1 


If M > 0, the camber will be read from the geometry 
dataset in the following manner. 

M * 1 The camber values only (not the twist or angle 

of attack) are to be read from the geometry 
dataset. Any twist in the initial camber is 
removed. 

M * 2 The camber values only (not the twist or angle 

of attack) are to be read from the geometry 
dataset. 

M ■ 3 The twist and camber values (no angle of attack) 

are to be read from the geometry dataset. 

M * 4 The twist, camber and angle of attack are to 

be read from the geometry dataset. 


If F i 0, the camber will be read from the geometry 
dataset in the following manner. 

F ^ 2 The flap values will be read from the geometry 

dataset and subtracted from the initial camber. 

F * 4 After subtracting from the initial camber, any 

flap deflections from the geometry dataset 
are set equal to zero. 
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method # 2 

Unless N < 0, twist and camber values may be obtained 
from normal velocities (normal facing upward from the 
panel) input in the order of the control points 
beginning in locations 2840 and 3340. The values input 
in these locations may include both twist and camber. 

If no values are input, there will be no be no camber 
input by this method, since all data locations are 
initially set equal to zero for the first case. The 
values from the 2840's are assumed to represent a 
smooth camber distribution and will be interpolated to 
compute drag. The values from the 3340 's are assumed 
to be due to flap deflections and will not be 
interpolated when drag is computed (see DATA(24) ) . 


method # 3 

If N > 0 camber may be input in the following manner. 

N ■ ND > 2, the number of x/c values where the z/c values 
for camber are specified (locations 3951,4001). These 
x/c stations are the same for all span stations and 
must include x/c * 0. and x/c » 1. The camber, and 
z/c's apply to fixed span stations only (see DATA(81) ) , 
unless DATA( 17) > 0. 


31 XOO The x value of the pivot point for angle of attack. 

Used for second order theory only. 

32 XCG The desired center of pressure (to override value from 

the geometry dataset. If a value < -99999. is input the 
XCG constraint is not used. 

33 ICG For non-symmetric rolling moment only (see DATA(ll)) 

34 ZCG For non-symmetric rolling moment only (see DATA(ll)) 

35 The freestream Mach number (used for computing leading 

edge suction only). 

If 0. the value from the geometry dataset will be used. 

If < 0. the Mach number will be set equal to 0. 

36 RC The x/c on each panel where the normal velocities (local 

angles of incidence) are interpolated for optimizing drag 
and for computing drag. 

If RC < 0. OPT @ x/c = - RC analysis @ x/c * 0.515 

If RC * 0. OPT § x/c • 0.875 analysis € x/c ■ 0.515 

If RC > 0. OPT § x/c * 0.875 analysis § x/c * RC 

Recall that the drag increment from each panel is equal 
to the delta Cp times the area times the local angle of 
incidence. The local angle of incidence may vary in 
the chordwise direction, and the boundary condition is 
satisfied on each panel only at x/c = 0.875. Therefore 
the interpolation is necessary to obtain an angle of 
incidence for each panel which is a closer 
approximation to the average value for the panel. 

37 RS The x/c fraction of each panel which is used to 

curvefit (x,Cp(x)) in order to obtain an approximate 
value for the leading edge suction (default * 0.250). 

Since RS < 1.0, if the value is input as M.RS, with 

M i 1, M will be the number of functions used for the 

curvefit. If M » 0, the default value is 3 functions 
(i.e. M = 3). If < 0. a detailed printout results. 

Drag polar (41 points given) 

38 Delta CL for drag polar. Default * 0.05 

39 Starting CL for drag polar. Default * 0.0 
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40 CBAR 


The reference chord length for computing Cm 
If 0. CBAR « CAVG is used. 

41 CAVG The reference chord length. 

If 0. value from the geometry dataset is used. 

42 SREF The reference area. 

If 0. value from the geometry dataset is used. 

43 Span Span, any consistent set of units may be used. This 

value is used for reference purposes only. 

If 0. value from the geometry dataset is used. 

50 Initial guess for angle of attack (for configurations 

having bodys paneled with source panels) 

51 ALPHAO(l) The configuration angle of attack with respect to the 

freestream (degrees). This is in addition to any 
camber distribution. For an optimization the 
resulting angle of attack is used for ALPHAO(l). If 
a value is input for ALPHAD(l) , and a body is 
present, it will be assumed to be a constraint on 
the body angle of attack. 


52 ALPHAD(2) The 2nd angle of attach where an analysis is desired. 
58 ALPHAD(8) (maximum number of angles of attack =8) 

If ALPHAD(K) < -90. ALPHAD(K-l) is used. 

If ALPHAD(K) > 90. an add load solution is given. 

i.e. (Alpha = 1.0) - (Alpha = 0.0) 

61 CAMTHK(l) Input in form OB-OV.CAM.THK, and used with ALPHAD(l) 

This input is used to define which set of unit 
solutions to combine into a final solution. Since 
the solution is linear, the unit solutions may 
be combined in any combination. 

OB > the order of the Cp on the body (one digit). 
OV = the order of the Cp on the lifting surfaces). 
CAN > 1 camber included. 

THK > 1 thickness included. 

CAM = 0 no camber included. 

THK * 0 no thickness included. 

OB, OV, CAM, THK are each one digit. 
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If=0. OB * 4, is used for the body Cp, 

OV, is determined by DATA(8) , and 

camber and thickness are included. 

OV = 2, if DATA(8) < 0. 

If < 0. Same asabove except the 2nd order axial 
solution (if performed) is used to 
calculate u,v,w. 

The body pressure formula is 
determined by the value of OB. 

OB Cp 

1 Cp = - 2 * u 

2 Cp * - 2 * u ~ beta2 * u*u - v*v - w*w 

3 Cp = - 2 * u - u*u ~ v*v - w*w 

4 Cp ■ Isentropic pressure formula (default) 

5 Cp ■ Isentropic for Alpha « 0. + isentropic add load. 

u,v f w Use 1st order axial contributions if camthk > 0 
u,v f w Use 2nd order axial contributions if camthk < 0 

E.g. CAMTHK = - 31.02 

Indicates, 2nd order u,v,w from axial solution (if performed). 

3 indicates, pressure formula f3 on body. 

1 indicates, 1st order Cp on lifting surfaces. 

0 indicates, camber is not included. 

2 indicates, thickness is included. 


62 CAMTHK(2) Used with ALPHAD(2) 
68 CAMTHK(8) Used with ALPHAD(8) 


The desired values of CL for each surface, (maximum of 10) 
If 0. the CL will not be constrained for that surface. 
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81 The input depends on whether a Cp optimization, a twist 
optimization, or an analysis is to be performed. 


For a Cp optimization ( DATA( 17) = 0. ) * * * 


The input consists of the list of span station numbers where 
a fixed camber is desired, followed by a list of span 
stations (input as negative numbers) where a fixed Cp 
distribution is desired. 

The span stations where a fixed camber is desired are 
referred to as fixed span stations, and the resulting 
camber slopes will be equal to the initial input camber. 

The list of span station numbers must be in ascending 
order. 

These span stations have the following properties: 

1. If a span station number in the list is input as J.R, 

R > 0, rather than just J., then the J'th span station 
is assumed to be a fixed span station with an unknown 
value of twist (see 351 for a definition of twist). 

All subsequent fixed span stations will have this value 
of twist (plus any twist difference specified by the 
inital input camber, which may contain twist, or the 
values obtained from locations beginning at 351), until 
another J.R is found in the list of fixed span 
stations, (e.g. this option could be used for a fixed 
surface with an unknown deflection). 

2. There are no Cp constraint functions used. Each panel Cp 

is regarded as an unknown (this can be thought of as a 
delta function type constraint function). 

3. The initial twist and camber is used (see the discussion 

for data location 29). 

4. The final slopes are obtained by adding the optimized 

angle of attack to any initial input twist and 
camber. 

5. No x/c-C.P. constraint is allowed. 

6. No CL*c/CAVG constraint is allowed. 
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The span stations where a fixed Cp distribution is desired 
are input as negative numbers following the list of span 
stations where a fixed camber is desired. The list of span 
station numbers must be in ascending order. The Cp 
distribution for this span station is input in place of 
normal velocities in locations beginning at 2840. 

i.e. If the span station has panel numbers beginning with 
ijl and ending with ij2, the Cp's are input in locations 
(2840 + ijl - 1) thru (2840 + ij2 - 1). 

These span stations have the following properties: 

1. There is one Cp constraint functions used at each of 

these stations. This constraint function has the 
desired Cp distribution. 

2. If a span station number in the list is input as - J.K, 

K > 0, rather than just - J., then the J'th span 
station will have the desired Cp distribution 
multiplied by a constant to be determined during the 
optimization. 

3. No x/c-C.P. constraint is allowed. 

4. No CL*c/CAVG constraint is allowed. 

5. Twist constraints may be used in the same manner as 

with span stations which do not have a fixed camber 
(see the above discussion of fixed span stations and 
the discussion of twist constraints at data location 
351). 


For a twist optimization ( DATA(17) > 0. ) * * * 

These locations are used to define the panels where flaps 
are located. It may be recalled that a twist 
optimization is only a term used to designate the type 
of optimization procedure involved. The twist 
optimization procedure may be used to find optimum 
twist camber and flap deflections. Further properties 
of these flaps (deflection ratios) may be input in 
locations beginning with 3340. 
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1. An input value of J.MMNNKK means that for span station 
J flap number KK runs from panels MM to panels nn. 

If KK is input as 00 it will be set equal to 01. 

A maximum of 20 flaps are allowed (KK=20). 


* * For an analysis ( DATA( 12) = 1. ) * * * 

For an analysis these locations are used to specify the 
span stations where the Cp's are to be set equal to 
zero, or the normal velocities on all panels of the 
span station are to be set equal to zero. 

This feature may be used to find the normal velocities 
induced on any set of span stations (or surfaces) by 
the remaining set of span stations (or surfaces). The 
normal velocities induced at the panel control points 
by all span stations where the panel Cp's were not set 
equal to zero may be found under the heading "DRAG 
N- VELOCITIES" in the printed output. Angle of attack 
and the presence or absence of camber may be controlled 
using locations 51-58 and 61-68. Data location 36 
should be set equal to 0.875 to supress the 
interpolation of the normal velocities to other than 
the panel control point (x/c = 0.875). The average 
incidence of induced velocities, for the first angle 
of attack only, may be found under the twist column 
in the printed output, or in the plot output dataset 
as twist. 

The span station numbers should be input in ascending 
order. 

1. A span station number input as a positive number will 

cause the panel Cp's to be set equal to 0.0 on that 
span station. The panel Cp's on all remaining span 
stations (except those as described below) will remain 
at the value achieved after an analysis at the 
specified angle of attack, and with or without the 
presence of camber. 

2. A span station number input as a negative number will 

result in the normal velocities on each panel of that 
span station to be equal to 0.0. The Cp's on all 
span stations input as negative numbers are adjusted 
to achieve this result by canceling the normal 
velocities induced by all panels where the Cp's were 
not set equal to 0.0. 
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151 


** For a Cp optimization (DATA(17) =0.) ** 

The desired values of CL*C/CAVG at each span station. 

If 0. is input, the value is unconstrained. 

** For a twist optimization (DATA(17) >0.) ** 

Used in conjunction with DATA(5), these locations are 
used to fix the coefficients of the camber constraint 
functions to desired values (input < 1.), or to constrain 
the value to be equal to the value at another span 
station (input the other span station number l 1.). if 
left 0., the coefficient value is unconstrained and an 
optimized value will be chosen. Only the camber 
constraints are input at these locations. For twist 
constraints see location 351. 

If a camber constraint coefficient of a given order is 
constrained, all higher order coefficients at the same 
station are also constrained, even if zero is input. 

Input is in order of the functions at each span station. 

151 1st camber function for station # 1 

152 2nd- camber function for station # 1 

150+(NFX-1) 1st camber function for station # 2 

2nd camber function for station # 2 


Locations 201 - 350 are used for Cp optimizations only. 

201 The relative value of delta CL at each span station 
If no values are input, the constraints will be 
satisfied exactly. 

251 The desired values of x/c-C.P. At each span station 
If 0. is input, the value is unconstrained 

301 The relative value of delta x/c-C.P. at span stations 

At least one of these values must be nonzero if any x/c-C.P. 
values are asked for. If a uniform change in x/c-C.P. is 
desired the value in these locations should be equal to the 
expected value of CL*c/CAVG. 
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351 These locations are used for inputing initial twist values or 
twist constraints. The input in these locations depends on 
whether a Cp optimization, a twist optimization, or an 
analysis is being performed. The twist of a span station is 
defined below. 


* * * For a Cp optimization ( DATA (17) = 0. ) * * * 


If span station J is a fixed span station, the values 
input in DATA(350+J) are treated as part of the initial 
twist input for span station J. Therefore if span 
station JP is the immediately preceeding fixed span 
station (to span station J), which had an unknown value 
of twist (J.K in DATA(81-150)), the resulting 
optimization will have: 

TWIST(J) = TWIST(JP) + ( TWISTO(J) - TWISTO(JP) ) 

where TWISTO(J) and TWISTO(JP) are the initial twist 
input values for span stations J, and, JP 
respectively (see data(29) for the input of initial 
twist and camber. 

If J - 1 the angle of attack of the first span station 
will be equal to DATA(351) plus any initial angle of 
attack (from locations beginning at 2840 and 3340 and 
the geometry input dataset (fixed span stations only). 

If span station J is not a fixed span station, any values 
input in DATA(350+J) are will be used to constrain the 
twist difference between span station, J, and a 
reference span station, JR. If DATA(350+J) = 0. the 
value of TWIST(J) will be unconstrained. Therefore the 
resulting optimization will have: 

TWIST(J) = DATA ( 3 50+ J) + TWIST(JR) 

JR = DATA ( 950+ J) 

if JR ■ 0 JR = J-l is used 

if JR > 0 JR ■ JR is used 

if JR < 0 the constraint is made with respect 

to the free stream, i.e. 

TWIST(J) = DATA ( 350+ J) 
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* * * 


For a twist optimization ( DATA(17) > 0. ) 


* * * 


If DATA(350+J) * 0., then the twist of span station J will 
not be constrained. 

Any values input in DATA(350+J) will be used to constrain 
the twist difference between span station, J, and a 
reference span station, JR. If the initial twist of 
station J was set equal to zero (see data(401)), the 
twist constraint is made on the final twist difference 
between the stations. If the initial twist was not set 
equal to zero, the constraint is made on the difference 
in the twist increments which are added to the initial 
twists at the two span stations. 

final twist difference constraint 

TWIST(J) * DATA ( 350+ J) + TWIST(JR) + TWISTO(JR) 

constraint on increment 

TWIST(J) = DATA(350+J) + TWIST(JR) 


Where: 

TWISTO(L) is the initial input twist of any station L 
TWIST(L) is the increment of twist give to station L 

TWIST(L) + TWISTO(L) is the final twist of station L 


JR = DATA ( 950+ J) 

if JR = 0 JR * J-l is used 

if JR > 0 JR * JR is used 

if JR < 0 the constraint is made with respect 

to the free stream, i.e. 

TWIST(J) - DATA(350+J) 

If J * 1 the angle of attack of the first span station 
will be equal to DATA(351) plus any initial angle of 
attack (from locations beginning at 2840 and 3340 and 
the geometry input dataset. 
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* * * 


For an analysis 


( DATA(12) = 1. ) 


* * * 


TWIST(J) » DATA(350+J) + TwistO(J) 

i.e. the initial twist of span station J is 

Increased by a fixed increment ( DATA ( 350+ J)). 

The twist of a given span station is defined as the angle of 
attack of the span station (with any camber removed) with 
respect to the freestream, when the angle of attack of the 
configuration is zero. This means the twist of span station J 
is equal to the local angle of attack with respect to the 
freestream, A0(J), reduced by the angle of attack, Alpha, 
times the local dihedral angle, COSZ(J). 

i.e. TWIST(J) = A0(J) - Alpha * COSZ(J) 

Alpha Is either the angle of attack of the reference 
span station or zero (see DATA(21>). 

The angle of attack. Alpha, may be obtained from the local 
angle of attack, A0(K), of the reference span station K, 
which is measured with respect to the freestream. 

Alpha = - A0(K) / COSZ(K) 

e.g. Suppose reference span station K and span station J 
have local angles of attack A0(K), and A0(J), with 
both A0(J) and A0(K) measured with respect to the 
freestream. Let COSZ(K) and COSZ(J) be the local 
dihedral angles. Then if a(K) and A(J) are the 
respective angles of attack after a pitch angle Alpha: 

A(K) = A0(K) + Alpha * COSZ(K) 

A(J) = A0(J) + Alpha * COSZ(J) 

The twist of span station J is defined as the local angle 
of attack, A(J), when the reference angle of attack, 

A(K) = 0. Therefore using the above: 

Alpha = - A0(K) / COSZ(K) , and 

TWIST(J) = A0(J) - A0(K) * COSZ(J) / COS(K) 
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401 ** Fcr a Cp optimization ** 

The span stations where no constraint functions are 
desired. See also location 5. 

** For a twist optimization ** 

Input as J.KL, where J is the span station where the camber 
or twist from the geometry dataset are to be set equal to 
zero (see location 29 for a discussion of input camber). 

If K > 0 the camber values are set equal to zero. 

If L * 0 the twist values are set equal to zero. 

e.g. an input value of 12.02 means that span station 12 
will have the camber, but not the twist, from the input 
dataset set equal to zero. 

The J values (span stations) must input in ascending 
order. 


551 The x/c values where the camber is specified (see 601). 

601 The z/c values for camber. 

The values for x/c = 0. and 1. must be specified. 

The z/c values for each x/c at every span station are 
input in consecutive locations starting with location 
451. The values are input first in the order of the 
x/c values for a given span station, and next in the 
order of the span stations. Any location left blank 
will be assumed to have a z/c value of 0. 

Camber values should be input only for fixed span 
stations (location 81), or when a twist optimization 
is preformed (location 17). Twist values may be 
input using locations 351 to 400. 

If an analysis only is being performed, data location 17 
should be made > 0. In this case all span stations must 
be input (although 0. is permissible). This is to avoid 
confusion of the fixed span stations with span stations 
where the Cp is to be set equal to zero (see the 
discussion for an analysis under location 81). 
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951 Twist constraint stations. 

L = DATA ( 950+ J) 

Twist constraint at span station J implies a 
nonzero of DATA(350+J). 

L < 0 Any twist constraint at span station j will be 
made with respect to the free stream. 

L = 0 Any twist constraint at span station j will be 
made with respect to span station (j-1). 

If j*l the constraint will be with respect to 
the freest ream unless a body is present. 

L > 0 Any twist constraint at span station j will be 
made with respect to span station L. 

L > nst Any twist constraint at span station j will be 
made with respect to the body. 


Hinge Moments locations 1001 - 1100 

The values in the following locations are for calculating or 
constraining up to twenty-five hinge moment coefficients. The 
hinge moments are calculated by taking the sum of the moments 
induced about a given line by a set of prescribed panels. 

Hinge moments may be constrained only when performing a twist 
optimization (DATA(17) >0.). 

The two points which determine each line about which the moments 
are taken, are input in locations beginning at 1001. The 
panels determining each integration area are input in 
locations beginning at 1051. 


1001 The z value of the 1st point for hinge moment line number 1. 

1002 The y value of the 1st point for hinge moment line number 1. 

1003 The z value of the 1st point for hinge moment line number 1. 

1004 The x value of the 2nd point for hinge moment line number 1. 

1005 The y value of the 2nd point for hinge moment line number 1. 

1006 The z value of the 2nd point for hinge moment line number 1. 

1007 CBAR for hinge moment number 1 (if 0. uses CAVG). 

1008 The constrained value for this hinge moment (no constraint if 0) 

1009 if < 0. hinge moment array number 1 will not be printed. 

1010 Reference area for hinge moment number 1 (if 0. uses SREF). 
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1011 The z value 

1012 The y value 

1013 The z value 


of the 1st point 
of the 1st point 
of the 1st point 


for hinge moment 
for hinge moment 
for hinge moment 


line number 2. 
line number 2. 
line number 2. 


1014 The x 

1015 The y 

1016 The z 


value of the 
value of the 
value of the 


2nd point for 
2nd point for 
2nd point for 


hinge moment 
hinge moment 
hinge moment 


line number 2. 
line number 2. 
line number 2. 


1017 CBAR for hinge moment number 2 (if 0. uses CAVG). 

1018 The constrained value for this hinge moment (no constraint if 0) 

1019 if < 0. hinge moment array number 2 will not be printed. 

1010 Reference area for hinge moment number 1 (if 0. uses SREF). 


1051 Beginning in this location the panel numbers for determining 

each hinge moment are listed. An input value of 0. indicates 
the last panel number for a given hinge moment calculation has 
been input. A sequence of panel numbers may be input by 
inputing a negative panel number after a positive number. 

i.e. 1051 3 5 10 20 - 55 0 3-40 

will have panel numbers 3,5,10, and 20 thru 55 for #1 
and panel numbers 3 thru 40 for #2 

due to camber, input in same order as the control points. 
These values will added to camber values input by other 
means. These values are assumed to represent a smooth 
camber distribution and will be interpolated to comDute 
drag. 


1200 If > 0. The effect of the normal velocities due to thickness 

will be included in the optimization. 

1201 The normal velocities at the panel control points due to 

thickness. 

1801 The total weight distributed over each panel. These values 
are used in conjunction with location 10 to obtain the 
deflections induced by the inertial loads on a flexible 
configuration. 
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2840 The normal (outward) velocities at the control points 

due to camber, input in same order as the control points. 

These values will added to camber values input by other 
means. These values are assumed to represent a smooth 
camber distribution and will be interpolated to compute 
drag. 

These locations are also used to input a desired Cp distribution 
for a given span station (see the discussion for data(81)). 

i.e. If the span station with a desired Cp distribution has 
panel numbers beginning with ijl and ending with ij2, the Cp’s 
are input in locations (2840 + ijl - 1) thru (2840 + ij2 - 1). 


3340 The same as for 2840 except the values are assumed to be 

due to flap deflections and will not be interpolated when 
drag is computed, nor will the the deflections contribute 
to camber when computing twist. 

If flap panel ratios are desired for a twist and flap 
optimization, (DATA(17) >0.), the flap panel ratios are 
input in these locations. I.e. if panel IJ is a panel of 
flap # j, and it is desired that this panel will deflect 
0.60 units for each unit deflection of flap # j, then set: 

Data ( 3340+1 J-l) = 100.60 

The value of 100. is added to differentiate between initial 
deflections and flap panel ratios. It is assumed that 
panel IJ is defined as a flap panel through locations 81-160. 
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LIST OF VARIABLES 


OPTIMIZATION VARIABLES 


NFX 

NFJ 

ISFX 

ISPJ 


I TRIM 
JTRIM 
KTRIM 


NCLS 

NCLY 

NPB 

NSTB 

NTWIST 

NDTWST 

NXCP 

CAM ( I J ) 
WF(IJ) 

Ml 

M2 

nu 

MLT 


The number of constraint functions at each station 
The number of constraint functions at each station (jet) 

The surface number where the jet attaches 
(if ISFX * NSF the jet is removed) 

The surface number of the jet 
(if ISFJ .gt. NSF) There is no jet 

The character array number for trim variation variable 
The constraint equation number for the trim variation. 
The variable which determines the desired variable for 
trim variation. 


The number of 
The number of 
The number of 
The number of 
The number of 
The number of 
The number of 


surface CL constraints 

spanwise CL constraints 

panels on fixed span stations 

span stations without constraint functions 

twist constraints 

unknown twists on fixed span stations 
x/c-C.P. constraints 


The normal velocity at the control point due to camber 
The deflections at the control points due to flaps. 

9 

The number of equations to be solved exactly 
The number of equations to be solved least square 
The total number of unknowns 
The total number of B(ML) 


GENERAL VARIABLES 


CHORD(J) 

IJS(IS) 

IJO(J) 

ISO(ISF) 

JS(IS) 

MS(IS) 

NB 

NCHORD(J) 

NS 

NSF 

NSPAN(ISF) 

NSL 

NST 


The chord value at the centroid of span station J 
The value of IJ where section IS begins 
The number of the first panel at span station J 
The section number where surface ISF begins. 

The span station where section IS begins 
The number of span stations in section IS 
The number of basic solutions. 

The number of panels at span station J 

The total number of sections 

The total number of surfaces 

The number of span stations on surface ISF 

The number of vortex shell sections 

The total number of span stations on all lifting 
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surfaces. 

NSTS The number of source span stations 

NTB The number of body panels. 

NTL The number of vortex shell panels. 

NTP The total number of panels. 

NTS The number of source parameters. 

NTSL The total number of vortex shell span stations 

NTV The number of vortex panels. 


LIFTING SURFACE PANELS 


CAM(IJ) 
COSZ(IJ) 
CX(I,IS) 
DS(J) 
DWDX(IJ) 
DX(I , IS) 
ETA(J) 

PA(IJ) 
SINZ(IJ) 
TWIST(J) 
WF(IJ) 
WTHK(IJ) 
X(KL,IC) 
XC(IJ,1) 
XC(IJ,2) 
XLE(1,J) 
XLE(2, J) 
XM(l f IS) 
XTE( 1 , J) 
XTE(2, J) 
XX(I,IS) 

y(kl,ic) 

YC(IJ) 
YCC(J) 
Y0(1,J) 
Y0(2,J) 
Z(KL,IC) 
ZC(IJ) 
ZCC(J) 
ZTHK(IJ) 
Z0(1,J) 
Z0(2, J) 


The normal velocity at the control point due to 
The cosine of the normal of panel IJ 
x/c values for the control points of section IS 
Width of span station J 

The local dv/dx at thickness control point IJ 

Delta (x/c) for the panels of section IS 

Fraction of span running distance for (YCC(J),ZCC(J)) 

camber. 

The area of panel IJ 

The sine of the normal of panel IJ 

Tvist of span station J 

The deflections at the control points due to flaps. 
The local dz/dx at thickness control point IJ 
x values of the corners of vortex panel KL (4) 
x value of the centroid for vortex panel IJ 

x value of the control point for vortex panel IJ 

x value at the leading left edge of span station J. 

x value at the leading right edge of span station J. 

x/c values for the midpoints of panels of section IS 
x value at the trailing left edge of span station J. 

x value at the trailing right edge of span station J. 

x/c values for the panels of section IS 
y values of the corners of vortex panel KL (2) 
y value of the control point for vortex panel IJ 
y value at the control point of span station J 
y value at the left edge of span station J. 

y value at the right edge of span station J. 

z values of the corners of vortex panel KL (2) 
z value of the control point for vortex panel IJ 
z value at the control point of span station J 
The local z/c at thickness control point IJ 
z value at the left edge of span station J. 

z value at the right edge of span station J. 
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BODY PANELS 


B(I ) 

BT(IC,IJ) 

BETA2X 


IBB(R) 

IJB(K) 

NX(K) 

NY(K) 


xb(ic,kl) 

XBO(I) 
XINLET(IJ) 
XN( 1 ,IJ) 
XN(2,IJ) 
XN(3,IJ) 


XN(4) 

XN(5) 

XN(6) 


SQRT(DY(l)**2 + BBTA2*DX(I )**2) 1-1,4 

tans**2 + beta2 for each side of body panel IJ 

* !• velocity boundary conditions. 

* I.- Mach**2 mass flux boundary conditions. 

The body station number of the first station in body 
section K. 

The panel number of the first body panel in section K. 
The number of x body stations for body section K 
The number of panels at each body station (half) for 
body section K. 

x values of the corners of body panel KL (4) 

The x value of the center of body station I. 

mass flow coefficient of body panel IJ; 0. = impermeable 

x - component of body panel IJ normal 

y - component of body panel IJ normal 

z - component of body panel IJ normal 

XN( 1) / XN(8) 

XN(2) / SQRT(XN(2)**2 + XN(3)**2) 

XN(3) / SQRT(XN(2)**2 + XN(3)**2) 


XN(7) 
XN(8) 
XN(9) 
XP(IC,KL) 
X0(1,IJ) 
X0 ( 2 , I J ) 
X0(3,IJ) 
YB(IC,RL) 
YP(IC,KL) 
ZB(IC,KL) 

zp(ic,kl) 


SQRT(XN(2)**2 + XN(3)**2) / XN(8) 
SQRT(beta2*XN(l)**2 + XN(2)**2 + XN(3)**2) 
SQRT( XN(1)**2 + XN(2)**2 + XN(3)**2) = 

x values of the corners of body panel KL (4) 

x value of the centroid of body panel IJ 

y value of the centroid of body panel IJ 

z value of the centroid of body panel IJ 

x values of the corners of body panel KL (4) 

x values of the corners of body panel KL (4) 

x values of the corners of body panel KL (4) 

x values of the corners of body panel KL (4) 


2. * area 
planar 


planar 

planar 
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COMPUTED VARIABLES 


Cp(IJ,K) The delta-Cp across panel IJ from basic solution K. 

UK(IJ,K) The control point upper surface x velocity. 

VR(IJ,K) The control point upper surface binormal velocity. 
HK(IJ,K) The control point upper surface normal velocity. 

PK(IJ,K) The control point upper surface velocity potential. 

US(IJ) The thickness induced upper surface z velocity. 

VS(IJ) The thickness induced upper surface binormal velocity. 

WS(IJ) The thickness induced upper surface normal velocity. 

PS(IJ) The thickness induced upper surface velocity potential. 

The following variables are required for 2nd order solutions only. 

WSX(IJ) d/dz of WS(IJ) source normal velocity. 

WCX(IJ) d/dz of WK(IJ,2) camber normal velocity. 

RKX The number of 2nd order boundary condition solutions. 

>6 if there is no body. 

* 9 if there is a body. 

UBE(IJ'K) The even symmetry u velocities due to 2nd order b.c. 
UB0(IJ,K) The odd symmetry u velocities due to 2nd order b.c. 
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SUBPROGRAMS 

A * signifies the subroutine is from program WBODT 


SUBROUTINE 


FUNCTION 


MAIN 

ABMAIN 

ACMAIN 
AER02 * 

AER02B * 

AER02P * 

AER02V * 

AMINMX * 

BAINFL * 


B0D7GM * 
B007MX 

B0D7RD * 
B0D7W * 

CAMBER 
CLCM * 

CLOPT 

CNSTRn 

CP2 

DECRD1 * 
DISPL7 * 

DMXMVE * 
DSSPL7 * 

ENDREC * 


Sets the main array size for the program. 

Reads input data, computes geometry, and sets array 
sizes. 

Controls the flow of the program. 

Controls the program flow for the computation of 1st 
or 2nd order pressures and loads. 

Computes the u velocities due to 2nd order boundary 
Conditions. 

Computes the 1st or 2nd order pressures from the 
induced velocities. 

Computes the odd and even symmetry velocities induced 
.by the basic solutions. 

Finds the largest, smallest or largest absolute value 
of the elements of an array. 

Computes normal velocities after the Cp’s or normal 
velocities on specified span stations are set equal to 
zero. 

Computes and prints body geometry from input data. 

Performs calculations for optimizations with paneled 
bodies. 

Reads body panel geometry from APAS dataset. 

Finds the intersection of the body and the aerodynmaic 
surfaces. 

Computes camber from input. 

Computes the lift drag and moment characteristics of 
aerodynamic surfaces. 

Controls flow of calculation for a Cp optimization. 

Calculates the Cp constraint functions for a Cp 
optimization. 

Adds second order terms to the Cp’s for an twist 
optimization using second order pressures. 

Reads the input data. 

Prints arrays of characteristics at the panel control 
points of the aerodynmaic surfaces. 

Moves the elements of a double precision array. 

Prints arrays of characteristics at the panel corner 
points of the aerodynamic surfaces. 

Used with errset to check for APAS type influence 
matrix. 
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FIELD * 


FXDX3 * 

FLXSOL 

FXROLL 

GEOM * 

GEOMR 
GAUSS 
HSHLDR * 

IHTRP3 * 

INTRP4 * 

INTRPX * 

INTRPY * 

MATRXF * 
MATRXT * 
HAXCHK 
MTXADD * 
MTXMLT * 
MTXMVE * 
NOTZRO * 
ORTHO * 


PHIXY * 

PLOT * 

POLAR * 

PTNFRM 

QIJT 

QUADMX 

REGION 
SQFIT * 


Computes first order field point properties using 
previously calculated influence coefficients and 
the first order solution. 

Integrates data using a third order curve fit through 
the nearest four points. 

Directs the calculation for the add load and basic load 
and obtains the deflections due to inertia forces. 

Prints rolling moment coefficients for rolling moment 
flap optimizations on flexible configurations. 

Computes the geometric characteristics of the 
aerodynamic surfaces and panels. 

Reads input geometry from geometry dataset. 

Solves simultaneous equations using Gaussian elimination. 

Solves sets of linear simultaneous equations in a 
least square sense. 

Interpolates or differentiates data using a third order 
curve fit through the nearest four points. 

Interpolates or differentiates data using a third or 
fourth order curve fit through the nearest four points. 

Interpolates or differentiates properties chordwise on an 
aerodynamic surface using subroutine INTRP4. 

Interpolates or differentiates properties spanwise on an 
aerodynamic surface using subroutine INTRP4. 

Displays arrays of data. 

Displays arrays of data. 

Checks for variable which exceed maximum value. 

Adds multiples of two arrays. 

Multiplies two arrays. 

Moves the elements of an array. 

Checks to see if an array has any nonzero elements. 

Solves sets of linear simultaneous equations using the 
method of succesive orthogonalizations. A 
quasi-inverse matrix may be computed or, if 
previously computed, used to perform the solution. 

Calculates v (binormal) velocities on aerodynamic 
surfaces from phi (for 2nd order solution). 

Writes geometry and aerodynamic data on a disk unit 
for computer graphics. 

Calculates first order drag polar and aerodynamic 
coff icients. 

Calculates Cp's from the constraint function coefficients 
in a Cp optimization. 

Computes drag coefficients for thickness drag in an 
optimization using 2nd order Cp's. 

Uses the method of Lagrange multipliers to maximize a 
quadratic form subject to a set of constraint equations. 

Calculates hinge moment coefficients. 

Calculates leading edge suction using a least square 
Curve fit of Cp to cot(x/c). 
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TRIM 

TRANSQ 

TWOPT 

TWOPTP 

TWOPTQ 

TWST 

VTXDRG * 
VTXDR2 * 
VTXMAX * 
VTXLFT * 
WCNSTR 

WINTRP * 


Used for trim polar calculations. 

Performs the transformation used for spanwise constraints 
on twist optimization constraint functions. 

Controls flow of calculation for a twist optimization. 

Calculates panel Cp's which result from a twist 
optimization. 

Calls QUADMX and checks for optimization of quadratic 
form. 

Separates twist and camber from normal velocities at 
panel control points and computes camber z/c's. 

Calculates vortex drag in the trefftz plane. 

Calculates coefficients for vtxdrg. 

Subroutine used for optimization with vortex lift. 

Calculates vortex lift from leading edge suction. 

Calculates the camber constraint functions for a twist 
optimization (NFX >1). 

Interpolates velocities to various points on panels. 


file 


use 


5 

6 
8 
9 

10 

11 

12 

13 

14 


input data 
Print output 
scratch file 
scratch file 
scratch file 

storage of influence matricies and quasi~inverse 

plot output file 

APAS geometry input file 

file for output of trim variation calculation 


143 



FLOW DIAGRAM 


********** 

* AAMAIN * 

********** 


Program OPT 


********** 

* ABMAIN *<==> 
********** 


********** 

<==>* DECRD1 * 

********** 

f 

********** 

<==>* ARRAY * 

********** 


********** ********** 

* ACMAIN *<==>)<==>* GEOMR * 

********** ********** 

V 

********** 

<==>* GEOM * 

********** 


********** 

<==>* CAMBER * 

********** 




********** ********** 

<==>* FLXSOL *<==>* WBSOL *<= 


********** 

********** 

<==>* CNSTRN * 

********** 

r 

********** 

<==>* BODYMX * 

********** 


********** 


(Cp opt only) 


*********** 

=>* ORTHO * 

*********** 


144 


A 

N 

A 

L 

T 

S 

I 

S 


V 


V 


V 


TWIST OPTIMIZATION 
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TEST CASS 


Results for the aspect ratio 2.5, sixty-three degree leading edge 
sweep trapezoidal wing of figure 4 are presented in this section. The 
10 X 10 uniform chordwise and spanwise aerodynamic paneling used for 
optimization is indicated and is typical of planform graphics output of 
the program. 

A M = 2.0, CL * 0.10 trimmed second order minimum drag due to lift 
case of figure 5 is presented in the remainder of this section. A root 
chord twist constraint was imposed to remove the small disturbance apex 
singularity common to subsonic leading edge problems. This result is 
compared to first and second order optima as a function of the 
longitudinal stability parameter dCm/dCl or alternatively the pitching 
moment at zero lift, Cmo. The supersonic nonlinear small disturbance 
minimum drag results are the first published to the knowledge of the 
author. The best second order result for the prestent problem is 6% 
lower than first order optimization and occurs at approximately twice 
the stability level. The first and second order lifting efficiency of a 
flat plate of the same planform is shown for comparison purposes. The 
impact of twist and camber for the subsonic leading edge case under 
consideration is substantial. Additional results are presented in 
reference 3. 

Test case input is presented on pages 151 and 152. Detailed 
program output for this case is presented on pages 153 through 202 
Typical aerodynamic data graphics output is presented on pages 203 
through 216. 
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Figure 4. Three Dimensional Second Order Optimum 
Model Problem 


149 


ORIGINAL PAS.Z 2 
OF POOR QUALITY 



-dCM/dC L 


Figure 5 . First and Second Order Optima for a 
Subsonic Edge Condition 
M = 2.0, C L = 0.1 
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TEST CASE INPUT 


0F poor c 


C63-45 

OPT 

T/C = 

0.050 

AT 

M = 

4.0 

1ST 

NO CONSTRAINT 

XCG = 

FREE 

C63-45 

THICK OPT 

T/C = 

0.050 

AT 

M = 

4.0 

2ND 

TWIST, NF=3.3 

XCG = 

12.0 

63-45 

OPT 

T/C = 

0.050 

AT 

M = 

4.0 

2ND 

TWIST, NF=3.3 

XCG = 

12.0 


FREE FORMAT 

C < 0. IF NO JETFLAP 

2 - 1 . 

C CL-TOTAL 

4 0.1000 

C ** # CHORD-CONST F'S < 0. • NONE 

5 3.3 

C MATRICES FROM WBODY , INCLUDING INVERSE (-3.0 INCLUDES THICKNESS) 

7- 3.5 

C 2ND ORDER ANALYSIS < 0. 

C* 2ND ORDER OPT =,< -1. ; < -2. FOR THKDRG 

8 - 1.00 

C ADD LOAD OPT, =0. IF < - 2., GEOMETRY ONLY WILL BE PRINTED. 

C 12 1.0 

12 0.0 

C IF > 0., THE NUMBER OF CONSTRAINT EQUATIONS PRINTED. 

14-100.100 

C CAMBER PRINTOUT =3.0 

15 3.0 

C GEOMETRY PRINTOUT =2.0 

16 0. 

C ** TWIST OPTIMIZATION IF > 0. 

17 1. 

C MOST PRINTOUT = 2.0, NO UPPER LOWER CP < 0. 

19-1. 

C > 0. FOR A PLOT DATASET 

20 2.0 

C* EXTEND FIRST SURFACE TO CENTERLINE 

C 26 1.2 

C CAMBER SPECIFICATION ( < 0. FOR NO CAMBER INPUT ) 

29-1. 

C* XCG 

32 12.000 
NO XCG CONSTRAINT 
32 -99999. 

MACH NUMBER ' 

35 4.00 

C CONTROL PT FOR DRAG CALCULATION U OPTIMIZATION IF < 0.) 

36 0.875 

C CONTROL PT FOR SUCTION CALCULATION < 0 FOR CURVE FIT PRINTOUT 

37 0.250 
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ORIGINAL PAGE i§ 
OF POOR QUALITY 


C 


C 

C 


c 


CAVG SREF SPAN 

41 0.0 160.0 0.0 

ODD AND EVEN SYMMETRY U,V,W 

45 334 . 

ODD AND EVEN SYMMETRY D/DX U,V,W 

46 334 . 

ANGLE OF ATTACK 

51 0.0 - 99.00 99.00 

CAM-THK 

61 11.22 12.22 12.02 

FIXED SPAN STATIONS 

81 1 . 2 . 

CAMBER CONSTRAINTS 
151 1 . E — 08 
153 0.0 0.0 

155 0.0 0.0 


TWIST 

351 0.0 - 1.50 0.0 

951 0.0 1.0 1.0 


1 0 . 

1 - 1 . 


0.0 

12.02 


0.0 

0.0 


0.0 

11.02 


0.0 

0.0 
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